PENDUGAAN PARAMETER DERET WAKTU HIDDEN MARKOV HAMILTON *
BERLIAN SETIAWATY, YANA ADHARINI DAN HIRASAWA
Departemen Matematika Fakultas Matematika dan Ilmu Pengetahuan Alam Institut Pertanian Bogor Jl Meranti, Kampus IPB Darmaga, Bogor 16680 Indonesia ABSTRAK. Pendugaan parameter untuk model deret waktu Hidden Markov Hamilton (1994) dilakukan mengunakan Metode Maximum Likelihood dan pendugaan ulang menggunakan metode Expectation Maximization. Dari kajian ini diperoleh algoritma untuk menduga parameter model. Kata kunci: Rantai Markov, Hidden Markov, Deret waktu Hidden Markov, Metode Expectation Maximization.
1. PENDAHULUAN Tulisan ini merupakan kajian tentang pendugaan parameter untuk deret waktu Hidden Markov, khususnya model Hamilton (1994). Pendugaan parameter menggunakan metode Maximum Likelihood dan pendugaan ulangnya menggunakan metode Expectatition Maximization (Metode EM). Dari kedua metode tersebut kemudian diturunkan suatu algoritma yang dapat dipakai secara umum untuk menduga parameter model deret waktu Hidden Markov Hamilton (1994). Tulisan ini dimulai dengan definisi model deret waktu Hidden Markov beserta sifat-sifatnya. Pada bagian 3 dibahas Pendugaan Parameter model dan terakhir pada bagian 4 dibahas pendugaan ulang parameter dan algoritmanya.
*Tulisan ini merupakan bagian dari hasil penelitian yang didanai oleh Hibah Penelitian PHK A2 Departemen Matematika IPB tahun 2005 dan 2006
13
14
BERLIAN SETIAWATY, YANA ADHARINI DAN HIRASAWA
2. MODEL DERET WAKTU HIDDEN MARKOV Pasangan proses stokastik {( X k , Yk ) : k } yang terdefinisi pada ruang peluang (, F , P) dan mempunyai nilai pada S Y disebut model hidden Markov apabila {X k } adalah rantai Markov dengan state berhingga dan diasumsikan bahwa rantai Markov {X k } tidak diamati. Sehingga {X k } tersembunyi (hidden) di balik proses observasi Yk . Banyaknya elemen dari S disebut ukuran (orde) dari model hidden Markov.
Pada bagian ini dibahas model hidden Markov Hamilton (1994) yang merupakan deret waktu yang memenuhi persamaan
Yt c( X t ) ( X t )Yt 1 t
(2.1)
di mana: X t adalah rantai Markov dengan ruang state S 1, 2 dan p P 11 p 21
p12 p 22
merupakan matriks peluang transisinya, di mana
pij P{ X t j | X t 1 i} .
Yt adalah proses yang diamati dan bernilai skalar.
{ t } adalah barisan peubah acak yang saling bebas dan berdistribusi normal N(0, 2 ). c (c1,c2 ) dan (1 , 2 ) 2 , dengan konstanta real. c( X t ) cX t dan ( X t ) X t .
Jadi
model
dicirikan
oleh
parameter
c1 , c2 , 1 ,
dan
2
adalah
c1 , c2 , 1 , 2 , 2 , P . Dengan
menggunakan metode EM akan diduga parameter c1 , c2 , 1 , 2 , 2 , P dari data Y. Karena t ~ N 0, 2 bebas stokastik identik maka dapat diperoleh fungsi sebaran bagi t sebagai berikut. yt
F t yt P t yt 0 yt
0
t 0 2 1 exp dt 2 2 2
t 2 1 exp d t 2 2 2
(2.2)
JMA, VOL. 5, NO. 1, JULI, 2006, 13-26
15
Berdasarkan persamaan (2.2) diturunkan fungsi sebaran bagi Yt dengan cara sebagai berikut. FYt yt P Yt yt P c X t X t yt 1 t yt
P t yt cX t X t yt 1
yt c X t X t yt 1
0
t 2 1 exp d t . 2 2 2
Misalkan v yt c X t X t yt 1 , maka v
FYt yt 0
t 2 1 exp d t 2 2 2
dan v 1 fYt yt FYt yt exp 2 yt 2 2
2
y c y 1 t Xt X t t 1 exp 2 2 2
y c y 1 t Xt X t t 1 exp 2 2 2
v yt 2
1
2
(2.3) Misalkan Yt adalah medan- yang lengkap dan dibangun oleh Y1 , Y2 , Y3 , Yt . Karena X t merupakan rantai Markov 2 state maka terdapat 2 fungsi kerapatan peluang bagi Yt . Kumpulan fungsi kerapatan tersebut dalam vektor 21 dilambangkan dengan t . Sehingga diperoleh: f yt X t 1, Yt 1 ; f yt X t 2, Yt 1 ;
t
1 yt c1 1 yt 1 2 exp 2 2 2 . yt c2 2 yt 1 2 1 exp 2 2 2
(2.4)
16
BERLIAN SETIAWATY, YANA ADHARINI DAN HIRASAWA
Misalkan ˆt t 1 melambangkan vektor 21 di mana elemen ke-j pada vektor merepresentasikan P X t j Yt 1 ; dan
melambangkan perkalian elemen
per elemen vektor, maka P X t 1 Yt 1 ; f yt X t 1, Yt 1 ; t t 1 t P X t 2 Yt 1 ; f yt X t 1, Yt 1 ; P X t 1 Yt 1 ; f yt X t 1, Yt 1 ; P X t 2 Yt 1 ; f yt X t 2, Yt 1 ;
(2.5)
P X t 1 Yt 1 ; f yt X t 1, Yt 1 ; . P X t 2 Yt 1 ; f yt X t 1, Yt 1 ;
Berdasarkan persamaan (2.5) maka dapat dituliskan: P X t j Yt 1 ; f yt X t j , Yt 1 ;
P X t j , Yt 1 ; P yt , X t j , Yt 1 ; P Yt 1 ; P X t j , Yt 1 ; P yt , X t j , Yt 1 ; P Yt 1 ;
(2.6)
P yt , X t j Yt 1;
sehingga diperoleh 2
f yt Yt 1 ; P yt , X t j Yt 1 ; j 1 2
P X t j Yt 1 ; f yt X t j , Yt 1 ; j 1
P X t 1 Yt 1 ; f yt X t 1, Yt 1 ; P X t 2 Yt 1 ; f yt X t 2, Yt 1 ;
1 ˆt t 1 t , di mana 1 1 1 . Jika persamaan (2.6) dibagi dengan persamaan (2.7) maka diperoleh
(2.7)
JMA, VOL. 5, NO. 1, JULI, 2006, 13-26
P yt , X t j Yt 1 ; f yt Yt 1 ;
17
P yt , X t j , Yt 1 ; P Yt 1 ; P Yt 1 ; P yt , Yt 1 ;
P yt , X t j , Yt 1 ; P X t j , yt , Yt 1; P yt , Yt 1 ; P yt , Yt 1;
P X t j yt , Yt 1 ; P X t j Yt ; .
(2.8) Sehingga berdasarkan persamaan (2.6), (2.7) dan (2.8) diperoleh
P X t j Yt ;
ˆt t
ˆ 1 ˆ
t t 1
P yt , X t j Yt 1 ;
t
t t 1
t
f yt Yt 1 ;
.
(2.9)
ti1 t P X t 1 i Yt ; 2
P X t 1 i, X t j Yt ; j 1 2
P X t 1 i X t j , Yt ; P X t j Yt ; j 1 2
P X t 1 i St j; P X t j Yt ; j 1 2
pijtt i
j 1
dan
t 1 t
p 1 p 2 11 t t 12 tt p11 1 2 p21 p p 21 22 tt tt
1 p12 t t P t t . p22 2 t t
(2.10)
Fungsi log likelihood L untuk data yang diamati Yt dapat dihitung dengan cara
L log f y1 Y 0 ; f y2 Y1 ;
f yT YT 1 ; log f yt Y t 1 ; T
t 1
(2.11)
18
BERLIAN SETIAWATY, YANA ADHARINI DAN HIRASAWA
Salah satu pendekatan yang dapat digunakan untuk memilih nilai awal bagi ˆt t 1 adalah dengan membuat 1 0 sama dengan vektor peluang tak bersyarat
1 2 yang memenuhi sifat ergodic, yaitu P 1 2 1 Sehingga diperoleh: T
1 p22 1 P X t 1 Y0 ; 2 p11 p22 ˆ 1 0 . 2 P X t 2 Y0 ; 1 p11 2 p p 11 22
Berdasarkan persamaan (2.7) diperoleh
f yt Yt 1 ; P X t 1 Yt 1 ; f yt X t 1, Yt 1 ; P X t 2 Yt 1 ; f yt X t 2, Yt 1 ; P X t 1 Yt 1 ; P X t 2 Yt 1 ;
xt c1 1 yt 1 2 1 exp 2 2 2 xt c2 2 yt 1 2 1 exp . 2 2 2 (2.12)
3. PENDUGAAN PARAMETER Penduga kemungkinan maksimum bagi diperoleh dengan memaksimum-kan fungsi log-likelihood T
L log f yt | Yt 1; . t 1
Berdasarkan persaman (2.12) diperoleh f yt Yt 1 ; c1
P X t 1 Yt 1;
yt c1 1 yt 1 2 yt c1 1 yt 1 1 exp 2 2 2 2
P X t 1 Yt 1; f yt X t 1, Yt 1;
yt c1 1 yt 1 . 2
Untuk memperoleh nilai cˆ1 yang memaksimum fungsi log-likelihood maka turunan pertama dari L terhadap c1 harus sama dengan nol.
JMA, VOL. 5, NO. 1, JULI, 2006, 13-26
19
L 0 c1 T
t 1
f yt Yt 1 ; 1 0 c1 f yt Y t 1 ;
T P X 1 Y ; f y X 1, Y ; t t t yt c1 1 y t 1 0 t 1 t 1 2 f yt Yt 1 ; t 1
T
P X t 1 Yt 1; f yt X t 1, Yt 1; f yt Yt 1;
t 1
T
c1
yt 1 yt 1
P X t 1 Yt 1; f yt X t 1, Yt 1; f yt Yt 1;
t 1
.
Sehingga diperoleh
P X t 1 Yt 1 ;ˆ f yt X t 1, Yt 1 ;ˆ yt ˆ1 yt 1
t 1
f yt Yt 1 ;ˆ
cˆ1
T
T
P X t 1 Yt 1 ;ˆ f yt X t 1, Yt 1 ;ˆ
t 1
f yt Yt 1 ;ˆ
(3.1)
(3.2)
Dengan cara yang sama diperoleh:
P X t 2 Yt 1 ;ˆ f yt X t 2, Yt 1;ˆ yt ˆ2 yt 1
t 1
f yt Yt 1 ;ˆ
cˆ2
T
T
P X t 2 Yt 1;ˆ f yt X t 2, Yt 1;ˆ
t 1
f yt Yt 1 ;ˆ
Berdasarkan persaman (2.12) diperoleh
f yt Yt 1; 1
yt c1 1 yt 1 2 yt c1 1 yt 1 1 P X t 1 Yt 1; exp 2 2 2 2 P X t 1 Yt 1 ; f yt X t 1, Yt 1 ;
yt c1 1 yt 1 yt 1 . 2
yt 1
20
BERLIAN SETIAWATY, YANA ADHARINI DAN HIRASAWA
. ˆ Untuk memperoleh nilai 1 yang memaksimum fungsi log-likelihood maka turunan pertama dari L terhadap 1 harus sama dengan nol. f yt Yt 1 ; L T 1 0 1 1 t 1 f yt Yt 1 ; T P X 1 Y ; f y X 1, Y ; t t t yt c1 1 yt 1 t 1 t 1 2 f yt Yt 1 ; t 1
yt 1 0
memberikan T
P X t 1 Yt 1 ; f yt X t 1, Yt 1 ; yt 1 yt c1 f yt Yt 1;
t 1
T
1
P X t 1 Yt 1 ; f yt X t 1, Yt 1 ; yt 1 f yt Yt 1 ;
t 1
Sehingga diperoleh
ˆ1
P X t 1 Yt 1;ˆ f yt X t 1, Yt 1;ˆ yt 1 yt cˆ1
t 1
f yt Yt 1;ˆ
T
f y Y ;ˆ
2 P X t 1 Yt 1;ˆ f yt X t 1, Yt 1;ˆ yt 1
t 1
2
.
T
(3.3)
t 1
t
Dengan cara yang sama diperoleh
P X t 2 Yt 1;ˆ f yt X t 2, Yt 1;ˆ yt 1 yt cˆ2
t 1
f yt Yt 1;ˆ
ˆ2
T
T
t 1
f y Y ;ˆ
2 P X t 2 Yt 1;ˆ f yt X t 2, Yt 1;ˆ yt 1
t
Berdasarkan persaman (2.12) diperoleh
t 1
(3.4)
JMA, VOL. 5, NO. 1, JULI, 2006, 13-26
f yt Yt 1 ; 2
21
y c y 2 1 t j j t 1 P X t j Yt 1 ; 3 exp 2 2 2 2 j 1 2
y c y 2 y c y 2 1 t j j t 1 t j j t 1 P X t j Yt 1; exp 2 2 2 4 2 2 y c y 2 1 t j j t 1 1 yt c j j yt 1 P X t j Yt 1 ; exp 2 2 2 2 2 4 2 j 1 2
2 1 yt c j j yt 1 P X t j Yt 1; f yt X t j, Yt 1; 2 2 2 4 j 1 2
2 y c y 2 t j j t 1 . 4 2 2 4
2
P X t j Yt 1; f yt X t j , Yt 1; j 1
Untuk memperoleh nilai ˆ 2 yang memaksimum fungsi log-likelihood maka turunan pertama dari L terhadap 2 harus sama dengan nol L 0 2 T
t 1
f yt Yt 1 ; 1 0 2 f yt Yt 1 ;
T 2 P X j Y ; f y X j , Y ; t t t 2 x c y 2 0 t 1 t 1 t St St t 1 f yt Yt 1 ; t 1 j 1 memberikan T 2 P X j Y ; f y X j , Y ; t t 1 t t t 1 2
f yt Yt 1 ;
t 1 j 1
T
2
P X t j Yt 1 ; f yt X t j , Yt 1; yt cSt St yt 1 f yt Yt 1 ;
t 1 j 1
Sehingga diperoleh T
2
ˆ 2
P X t j Yt 1 ;ˆ f yt X t j , Yt 1 ;ˆ yt cˆSt ˆSt yt 1
t 1 j 1
T
2
t 1 j 1
f yt Yt 1 ;ˆ
P X t j Yt 1 ;ˆ f yt X t j , Yt 1 ;ˆ
f yt Yt 1 ;ˆ
2
.
2
.
(3.5)
Untuk menduga parameter P digunakan formula (3.6) berikut yang diperoleh dari Hamilton (1990).
22
BERLIAN SETIAWATY, YANA ADHARINI DAN HIRASAWA
P X T
pˆ ij
t 2
T
t
j , X t 1 i | YT ;ˆ
P X t 1 i | YT ;ˆ t 2
(3.6)
di mana menurut Kim (1994)
P X t j , X t 1 i | YT ;ˆ
P X P X
j | Y ;ˆ P X i | X j , Y ;ˆ j | Y ;ˆ P X i, s j | Y ;ˆ P X j | Y ;ˆ P X j | Y ;ˆ P X i | Y ;ˆ P X i | X P X j | Y ;ˆ P X t j | YT ;ˆ P X t -1 i | X t j , YT ;ˆ t
T
t -1
t
T
t -1
t
t
t
t
T
T
T
T
t -1
T
t
t
t -1
i
T
(3.7) dan
2
P X t 1 i | YT ;ˆ P X t j , X t 1 i | YT ;ˆ . j 1
4. PENDUGAAN ULANG PARAMETER DAN ALGORITMA Karena persamaan (3.1) sampai (3.7) tak linear, maka tidak mungkin untuk memperoleh ˆ secara analitik sebagai fungsi dari Y1 , Y2 , Y3 , , YT . Sehingga untuk mencari penduga kemungkinan maksimum digunakan algoritma iteratif dari metode EM.
4.1 Metode Expectation Maximization (Metode EM) Algoritma EM dikembangkan oleh Baum and Petrie (1966) dengan ide dasar sebagai berikut. Misalkan
P :
adalah koleksi ukuran peluang yang terdefinisi pada
ruang (, G) dan kontinu absolut terhadap P0 . Misalkan Y G . Definisikan fungsi likelihood untuk menentukan penduga parameter berdasarkan informasi Y sebagai
JMA, VOL. 5, NO. 1, JULI, 2006, 13-26
23
dP L( ) E0 Y dP0 dan penduga maksimum likelihood didefinisikan sebagai ˆ arg max L( ) .
Secara umum penduga maksimum likelihood ˆ sulit dihitung secara langsung. Algoritma EM memberikan suatu metode iteratif untuk mengaproksimasi ˆ , dengan prosedur sebagai berikut. Langkah 1: Set p 0 dan pilih ˆ0 . Langkah 2: [Langkah-E] dP Set ˆp dan hitung Q , E log Y . dP
Langkah 3: [Langkah-M] Tentukan ˆp 1 arg max Q , .
p p 1 Ulangi langkah 2 sampai kriteria berhenti dipenuhi.
Langkah 4:
Catatan:
1. Barisan ˆp : p 0 memberikan barisan L ˆp : p 0 yang tak turun. 2. Menurut ketaksamaan Jensen, Q ˆp 1 ,ˆp log L(ˆp 1 ) log L(ˆp ) .
3. Q , disebut pseudo-loglikelihood bersyarat.
4.2 Pendugaan ulang parameter menggunakan metode EM Dengan menggunakan meode EM diperoleh algoritma untuk menduga ulang parameter model. Algoritma tersebut sebagai berikut. Langkah 1: Tentukan banyaknya data T yang akan diamati serta tentukan juga nilai
y0 , y1 , y2 ,
, yT 1 dan matriks transisi p P 11 p21
p12 . p22
Beri nilai awal bagi ˆ yang dilambangkan dengan ˆ( m ) cˆ1 , cˆ2 , ˆ1 , ˆ2 , ˆ 2 .
24
BERLIAN SETIAWATY, YANA ADHARINI DAN HIRASAWA
Langkah 2: Cari fungsi kerapatan bersyarat bagi yT untuk setiap t 1, 2, , T dengan cara
f y X 1, Y ;ˆ m t t t 1 t ˆ m f yt X t 2, Yt 1 ;
ˆ ˆ 1 exp yt c1 1 yt 1 2ˆ 2 2ˆ y cˆ ˆ y 1 t 2 2 t 1 exp 2 2ˆ 2ˆ
2
. 2
Langkah 3: Penarikan kesimpulan optimal dan peramalan untuk setiap waktu t pada contoh dapat diperoleh melalui iterasi: 3.1.
Tentukan nilai awal bagi ˆt t 1 yang dilambangkan dengan ˆ1 0 .
3.2.
Berikan nilai awal i 1 .
3.3.
Untuk t i , cari nilai dari
m f yt Yt 1 ;ˆ 1 ˆt t 1 t
m P X t 1 Yt 1; ˆ
P X t 2 Yt 1; ˆ
y cˆ ˆ y 1 t 1 1 t 1 exp 2ˆ 2 2ˆ
m
3.4. Ulangi mulai dari langkah (3.3). Stop jika t T . Lanjutkan ke langkah 4. Langkah 4: Cari nilai dari:
2
y cˆ ˆ y 1 t 2 2 t 1 exp 2ˆ 2 2ˆ
P X 1 Y ;ˆ m ˆt t 1 t t t-1 ˆ , t t m 1 ˆ ˆ t t t 1 P X t 2 Yt-1 ; P X 1 Y ;ˆ m t 1 t ˆt 1 t P ˆt t , m ˆ P X t 1 2 Yt ; i i 1 .
, 2
JMA, VOL. 5, NO. 1, JULI, 2006, 13-26
T
y y f y | Y ;ˆ ˆ f y | X i, Y ;ˆ f y | Y ;ˆ
ˆ j t 1|t f yt | X t i, Yt 1 ;ˆ
cˆi
j
t 1|t
t 1
ˆi
t 1
i 1, 2.
,
t 1
t
t
t
t 1
y cˆ y f y | Y ;ˆ ˆ f y | X i, Y ;ˆ y f y | Y ;ˆ
ˆ j t 1|t f yt | X t i, Yt 1 ;ˆ
t
j
2
t
t 1
t
t
,
t 1
f yt | Yt 1 ;
t 1 j 1
2
f y | Y ;ˆ
cˆ X ˆX yt 1 t
ˆ j t 1|t f yt | X t j , Yt 1 ;ˆ
t 1 j 1
i 1, 2.
t 1
ˆ j t 1|t f yt | X t j , Yt 1 ;ˆ yt T
2
t 1|t
t 1
T
t 1
i
t 1
t
T
ˆ 2
t 1
t
T
T
i t 1
t
t 1
25
t
t
2
.
t 1
Langkah 5: Beri nama parameter yang dihasilkan pada langkah 4 dengan ˆ cˆ , cˆ , ˆ , ˆ , ˆ 2 , k 0,1, 2, , T 1 . k 1
1
2
1
2
Langkah 6: Cari P yang baru, yaitu:
ˆt(|Tj ) ˆt(|tj ) P ˆt(j1|)T ()ˆt(j1|)t di mana () merupakan operasi perkalian elemen per elemen vektor, P X t j , X t 1 i | Y
P X t j | Y P X t -1 i | Yt P X t i | X t -1 i P X t j | Yt
ˆ j t |T ˆi t 1|t 1 pij , ˆ j t |t
26
BERLIAN SETIAWATY, YANA ADHARINI DAN HIRASAWA
2
P X t 1 i | Yt ;ˆ P X t j , X t 1 i | Y ,
P X T
pˆ ij
t 2
t
j , X t 1 i | YT ;ˆ
P X T
t 2
j 1
t 1
i | Y ;ˆ T
T
ˆ j t |T ˆi t 1|t 1 pij
t 2
T
2
t 2 j 1
ˆ j t |t . ˆ j t |T ˆi t 1|t 1 pij ˆ j t |t
Langkah 7: Selama k T , ulangi mulai dari langkah 2. Gunakan parameter yang sudah dihasilkan untuk mencari nilai harapan bagi nilai y yang akan datang.
E Yt 1 | X t 1 j, Yt ; c j j yt .
DAFTAR PUSTAKA [1]. Baum,L.E. and Petrie, T. 1966. Statistical inference for probabilistic functions of finite state Markov chains. Annals of the Institute of Statistical Mathematics, 37:1554-1563. [2]. Hamilton, J. D. 1990. Analysis of Time Series Subject to Changes in Regime. Journal of Econometrics 45: 39 -70. [3]. Hamilton, J. D. 1994. Time Series Analysis. Princeton University Press, New Jersey. [4]. Kim, C. J. 1994. Dynamic linear models with Markov-Swiching. Journal of Econometrics 60: 1- 22.