LAMPIRAN
23
Lampiran 1. Skema MCMC untuk model ARCH(1) non-central Student-t
1. Inisialisasi a0, b0, 0, dan z0. 2. Pembangkitan
.
Distribusi posterior bersyarat untuk diberikan oleh
.
Dalam kasus ini,
bisa dibangkitkan secara langsung dari distribusi normal yaitu
.
dengan
3. Pembangkitan
dan
.
Distribusi posterior bersyarat untuk
diberikan oleh
p( zt ) f IG zt | , t g ( zt | a, b, , Rt 1 , Rt ) , 24
(1 b) R12 a , t 1, 1 2 a dengan , t 2 R ( a bRt21 ) 2 t , t 2,...,T , 2( a bRt21 ) 1 (1 b) 2 R1 exp , t 1, 1 az12 g ( zt | a, b, , Rt 1 , Rt ) Rt exp , t 2,...,T . 1 1 2 2 2 a bR z t 1 t
Dalam hal ini, posterior zt tidak mengikuti suatu distribusi tertentu, maka parameter zt dibangkitkan dengan menggunakan metode IC-MH dengan proposalnya yaitu
zt* ~ IG( , t ) dan rasio penerimaannya yaitu: r zt* , zt
g ( zt* | a, b, , Rt 1 , Rt ) . g ( zt | a, b, , Rt 1 , Rt )
4. Pembangkitan Distribusi posterior bersyarat untuk v diberikan oleh . Diambil logaritma distribusi posterior bersyarat untuk :
Nilai
dibangkitkan menggunakan metode IC-MH, dengan proposal untuk dan probabilitas penerimaannya yaitu
Dicari modus posterior dan
dari
, artinya bahwa
yaitu
. . Dalam kasus ini
ditentukan dengan menggunakan metode yang didasarkan pada tingkah laku
distribusi di sekitar modus. Lebih lanjut diperoleh bahwa derivatif pertama dan kedua dari
berturut-turut yaitu:
25
Kemudian itu diambil
dan
. Namun, dengan
bisa bernilai positif, karena .
26
Lampiran 2. Skema MCMC untuk model ARCH(1) Skewed Student-t 1.
Inisialisasi a0, b0, 0, dan z0.
2.
Membangkitkan nilai acak a | b, , k , Z, R . Distribusi posterior bersyarat untuk a diberikan oleh
Diambil logaritma distribusi posterior bersyarat untuk a:
Nilai a dibangkitkan menggunakan metode IC-MH, dengan proposal untuk a yaitu dan probabilitas penerimaannya yaitu Dicari modus posterior
dari
derivatif pertama dan kedua dari
, artinya bahwa berturut-turut yaitu
27
. . Diperoleh bahwa
Kemudian itu diambil 3.
dan
. Namun dengan
bisa bernilai positif, karena .
Membangkitkan nilai acak b | a, , k , Z, R Distribusi posterior bersyarat untuk a diberikan oleh
Diambil logaritma distribusi posterior bersyarat untuk b:
28
Nilai b dibangkitkan menggunakan metode IC-MH, dengan proposal untuk b yaitu dan probabilitas penerimaannya yaitu
Dicari modus posterior
dari
, artinya bahwa
derivatif pertama dan kedua dari
Kemudian itu diambil
.
. Diperoleh bahwa
berturut-turut yaitu
dan
. Namun dengan
bisa bernilai positif, karena .
29
4.
Membangkitkan nilai acak k | a, b, , Z, R Distribusi posterior bersyarat untuk k diberikan oleh
Diambil logaritma distribusi posterior bersyarat untuk k:
dengan
dan
. Dalam kasus ini,
langsung dari distribusi normal yaitu
bisa dibangkitkan secara
dimana:
dan 5.
Membangkitkan nilai-nilai acak zt | a, b, k , , R Distribusi posterior bersyarat untuk
dengan
Nilai
diberikan oleh
,
, t=1,...,T
dibangkitkan menggunakan metode IC-MH, dengan proposalnya yaitu dan rasio penerimaannya yaitu
30
.
6.
Membangkitkan nilai acak | k , a, b, Z, R Distribusi posterior bersyarat untuk v diberikan oleh
Diambil logaritma distribusi posterior bersyarat untuk :
Nilai
dibangkitkan menggunakan metode IC-MH, dengan proposal untuk dan probabilitas penerimaannya yaitu
7. Kembali ke langkah 2
31
yaitu .
Lampiran 3. Kode Matlab untuk Estimasi Model ARCH(1) noncentral Student-t 3.1 Kode Utama function hasil=archnct_mcmc(Rt,HP) % % % % % % % % % % % %
Tujuan: Mengestimasi parameter-parameter dalam model ARCH(1) R_t = sigma_t * zt^0.5 * (miu+xhi_t), xhi_t~N(0,1) sigma_t^2 = a + b*R_{t-1}^2 ------------------------------------------------------------------------Algoritma: MCMC ------------------------------------------------------------------------Masukan: Rt = Returns = 100*[(log(S_t)-log(S_{t-1}))- 1/T*sum((log(S_t)log(S_{t-1}))], dimana S_t = nilai kurs pada saat t HP = Hyperparameter = [lambda alpha_b beta_b miu_c var_c alpha_nu beta_nu], untuk prior a ~ exp(lambda), b ~ beta(alpha_b,beta_b), c ~ Normal(miu_c,var_c), dan nu ~ Gamma(alpha_nu, beta_nu)
% ------------------------------------------------------------------------% Keluaran: hasil.draws = draws % ----- I: Inisialisasi T = length(Rt); R2 = Rt.^2; % prior untuk a lamd = HP(1); % prior untuk b alp = HP(2); bet = HP(3); %prior untuk c muc=HP(4); varc=HP(5); % prior untuk nu alpnu = HP(6); betnu = HP(7); % nilai awal a = 0.1; b = 0.1; nu = 20; zt = 1./gamrnd(nu/2,2/nu,T,1); % banyaknya replikasi Nits = 15000; BI = 5000; N = Nits-BI; % alokasi penyimpanan sampel av = zeros(N,1); bv = zeros(N,1); cv = zeros(N,1); nuv = zeros(N,1); zv = zeros(T,1); vol = zeros(T,1); % ----- Algoritma MCMC. Step 1: Membangkitkan Rantai Markov tic for its = 1:Nits % --- pembangkitan nilai acak mu Vc = 1/(T+1/varc); Mc = Vc*(Rt(1)*sqrt((1-b)/(a*zt(1)))... +sum(Rt(2:T).*((a+b*R2(1:T-1)).*zt(2:T)).^(-0.5))+muc/varc); c = normrnd(Mc,sqrt(Vc),1,1);
32
% --- pembangkitan nilai acak z alpz = (nu+1)/2; betz1 = 0.5*((1-b)*Rt(1)^2+a*nu)/a; betz2 = (R2(2:T)+(a+b*R2(1:T-1)).*nu)./(2*(a+b*R2(1:T-1))); betz = [betz1; betz2]; % pembangkitan proposal zprop = 1./gamrnd(alpz,1./betz); % mengevaluasi probabilitas penerimaan gzold1 = (sqrt(1-b)*c*Rt(1)./(a*zt(1))^(-0.5)); gzold2 = (c*Rt(2:T)./(sqrt(a+b*R2(1:T... -1).*sqrt(zt(2:T))))); gzold = [gzold1; gzold2]; gznew1 gznew2 gznew
= (sqrt(1-b)*c*Rt(1)*(a*zprop(1))^(-0.5)); = (c*Rt(2:T)./(sqrt(a+b*R2(1:T... -1).*sqrt(zprop(2:T))))); = [gznew1; gznew2];
ratio = exp(gznew-gzold); % pembaharuan u = rand(T,1); minrat = min(1,ratio); ind = find(u<=minrat); zt(ind) = zprop(ind); % mencari rata-rata dan variansi untuk proposal bersyarat hpv = [alpnu betnu]; mv = bisection_nu(hpv,zt); % rata-rata d2v = 0.5*T/mv-T/4*psi(1,mv/2)-(alpnu-1)/mv^2; Dv = min(-0.0001,d2v); s2v = -1/Dv; %variansi % algoritma IC-MH % pembangkitan proposal nu* pr = truncnormrnd(1,mv,sqrt(s2v),2.1,40); % mengevaluasi probabilitas penerimaan log_pv = 0.5*pr*T*log(pr/2)-T*gammaln(pr/2)-pr/2 ... *sum(log(zt)+1./zt)+(alpnu-1)*log(pr)-betnu*pr; post_pr = exp(log_pv); log_pv post_o
= 0.5*nu*T*log(nu/2)-T*gammaln(nu/2)-nu/2 ... *sum(log(zt)+1./zt)+(alpnu-1)*log(nu)-betnu*nu; = exp(log_pv);
ratio = post_pr/post_o; ap = min(1,ratio); % pembangkitan variabel acak seragam u = rand(1); % pembaruan if u <= ap, nu = pr; end % --- pembangkitan sampel a % mencari rata-rata dan variansi untuk distribusi proposal hp1 = lamd; m_a = bisection_nct(Rt,R2,a,b,c,hp1,zt,'a'); % rata-rata d2a = 1/(2*m_a^2)+0.5*sum(1./(m_a+b*R2(1:T-1)).^2)... -0.5*(1-b)*R2(1)/(m_a^3*zt(1))-(0.25*c*Rt(1)*(1... -b)^2./((1-b)/m_a*zt(1))^(0.5)*m_a^(4)*zt(1)^2)... +c*Rt(1)*(1-b)./sqrt((1-b)./m_a*zt(1))*m_a^(3)*zt(1)... -sum(R2(2:T)./((m_a+b*R2(1:T-1)).^3.*zt(2:T)))... +0.75*sum(c*Rt(2:T)./(m_a+b*Rt(1:T... -1).^(2.5).*sqrt(zt(2:T)))); Da = min(-0.0001,d2a); s2_a = -1/Da; % variansi % algoritma IC-MH % pembagkitan proposal a*
33
pr = truncnormrnd(1,m_a,sqrt(s2_a),1e-4,1); % mengevaluasi probabilitas penerimaan log_pa = -0.5*log(pr)-0.5*sum(log(pr+b*R2(1:T-1)))... -0.5*(1-b)*R2(1)/(pr*zt(1))... +sqrt(1-b/pr*zt(1)*c*Rt(1))... -0.5*sum((R2(2:T)./((pr+b*R2(1:T-1)).*zt(2:T)))... +(c*Rt(2:T)./(sqrt(pr+b*R2(1:T-1)).*sqrt(zt(2:T)))))... -lamd*pr; % F(pr) post_pr = exp(log_pa); log_pa
post_o
= -0.5*log(a)-0.5*sum(log(a+b*R2(1:T-1)))... -0.5*(1-b)*R2(1)/(a*zt(1))-... +sqrt(1-b/a*zt(1)*c*Rt(1)) ... -0.5*sum((R2(2:T)./((a+b*R2(1:T-1)).*zt(2:T))) ... +(c*Rt(2:T)./(sqrt(a+b*R2(1:T-1)).*sqrt(zt(2:T)))))... -lamd*a; % F(a) = exp(log_pa);
ratio = post_pr/post_o; ap = min(1,ratio); % pembangkitan variabel acak seragam u = rand(1); % pembaruan if u <= ap, a = pr; end % --- pembangkitan nilai acak b % mencari rata-rata dan variansi untuk distribusi proposal mup = alp; Vp = bet; hp1 = [mup Vp]; m_b = bisection_nct(Rt,R2,a,b,c,hp1,zt,'b'); % rata-rata d2b = -0.5/(1-m_b)^2+0.5*sum((R2(1:T-1).^2)./((a+m_b*R2(1:T... -1)).^2))-0.25*c*Rt(1)/((1... -m_b/a*zt(1)).^(1.5)*a^2*zt(1)^2)... -sum(((R2(1:T-1).^2).*(R2(2:T)))... ./((a+m_b*R2(1:T-1)).^3.*zt(2:T)))... +0.75*sum(R2(1:T-1).^(2)*c.*Rt(2:T)./(a+m_b*Rt(1:T... -1).^(2.5).*sqrt(zt(2:T))))... -(alp-1)/(m_b^2)-(bet-1)/(1-m_b)^2; Db = min(-0.0001,d2b); s2_b = -1/Db; % variansi % algoritma IC-MH % pembangkitan proposal b* pr = truncnormrnd(1,m_b,sqrt(s2_b),0,1); % mengevaluasi probabilitas penerimaan log_pb = 0.5*log(1-pr)-0.5*sum(log(a+pr*R2(1:T-1)))... -0.5*(1-pr)*R2(1)/(a*zt(1))... +sqrt(1-pr/a*zt(1))*c*Rt(1)... -0.5*sum((R2(2:T)./(a+pr*R2(1:T-1).*zt(2:T)))... +c*Rt(2:T)./sqrt(a+pr*R2(1:T-1)).*sqrt(zt(2:T)))... +(alp-1)*log(pr)+(bet-1)*log(1-pr); % F(pr) post_pr = exp(log_pb); log_pb
post_o
= 0.5*log(1-b)-0.5*sum(log(a+b*R2(1:T-1)))... -0.5*(1-b)*R2(1)/(a*zt(1))... +sqrt(1-b/a*zt(1))*c*Rt(1)... -0.5*sum((R2(2:T)./(a+b*R2(1:T-1).*zt(2:T)))... +c*Rt(2:T)./sqrt(a+b*R2(1:T-1)).*sqrt(zt(2:T)))... +(alp-1)*log(b)+(bet-1)*log(1-b); %F(b) = exp(log_pb);
ratio = post_pr/post_o; ap = min(1,ratio); % pembangkitan variabel acak seragam u = rand(1); % pembaruan if u <= ap, b = pr; end
34
% pengestimasian sigma volt = [a/(1-b); a+b*R2(1:T-1)]; % simpan a, b, c, dan nu if its > BI av(its-BI,1) = a; bv(its-BI,1) = b; cv(its-BI,1) = c; nuv(its-BI,1) = nu; zv = ((its-BI-1)*zv+zt)/(its-BI); vol = ((its-BI-1)*vol+volt)/(its-BI); end end toc % ----draws = MP = SP =
Algoritma MCMC. Step 2: Menghitung rata-rata Monte Carlo [av bv cv nuv]; mean(draws); std(draws);
% ===== Integrated Autocorrelation Time (IACT) ============================ % Berapa banyak sampel yang harus dibangkitkan untuk mendapatkan sampel % yang independent (seberapa cepat konvergensi simulasi) resultsIAT = IACT(draws); IAT = [resultsIAT.iact]; % ===== Uji Konvergensi Geweke ============================================ idraw1 = round(.1*N); resultCV = momentg(draws(1:idraw1,:)); meansa = [resultCV.pmean]; nsea = [resultCV.nse1]; idraw2 resultCV meansb nseb
= = = =
round(.5*N)+1; momentg(draws(idraw2:N,:)); [resultCV.pmean]; [resultCV.nse1];
CD = (meansa - meansb)./sqrt(nsea+nseb); onetail = 1-normcdf(abs(CD),0,1); pV = 2*onetail; % ===== 95% Highest Posterior Density (HPD) Interval ====================== resultsHPD = HPD(draws,0.05); LB = [resultsHPD.LB]; UB = [resultsHPD.UB]; % % ===== Numerical Standard Error (NSE) ==================================== resultsNSE = NSE(draws); NSEd = [resultsNSE.nse]; %====================== Pengaturan Pencetakan Hasil ======================= %----- Statistik Parameter in.cnames = char('a','b','c','nu'); in.rnames = char('Parameter','Mean','SD','LB','UB','IACT','NSE','G-CD','pValue'); in.fmt = '%16.6f'; tmp = [MP; SP; LB; UB; IAT; NSEd; CD; pV]; fprintf(1,'Estimasi menggunakan MCMC dan Uji Diagnostik\n'); mprint(tmp,in); hasil.vol hasil.zv hasil.av
= vol; = zv; = av;
35
hasil.bv hasil.cv hasil.nuv hasil.draws
= = = =
bv; cv; nuv; draws;
%========== Calculate Marginal likelihood by Gelfand-Dey (1994) =========== theta = draws'; [npar,ndraw] = size(theta); mvn_mean = mean(theta,2); mvn_cov = cov(theta'); inv_mvn_cov = inv(mvn_cov); %----- POSTERIOR % domain truncate_val = zeros(1,ndraw); for i = 1:ndraw draw_element = theta(:,i); truncate_val(i) = (draw_element-mvn_mean)' * inv_mvn_cov... * (draw_element-mvn_mean); end critical_value = chi2inv(0.99,npar); truncate_index = truncate_val > critical_value; source_pdf_log = mvn_pdf_log(theta,mvn_mean(:,ones(ndraw,1)),mvn_cov) ... - log(0.99); source_pdf_log(truncate_index) = -inf; %----- PRIOR % a prior_pdf_a_log = log(lamd)-lamd*av'; % b prior_pdf_b_log = -betaln(alp,bet)+(alp-1).*log(bv')+(bet-1).*log(1-bv'); % c prior_pdf_c_log = -1/2*log(2*pi*varc) - 1/2*(cv'-muc).^2/varc; % nu prior_pdf_nu_log = -gammaln(alpnu)+alpnu.*log(betnu)+(alpnu-1).*log(nuv')betnu.*nuv'; % LOG TOTALLY PRIOR prior_pdf_log = sum(prior_pdf_a_log) + prior_pdf_b_log ... + prior_pdf_c_log + prior_pdf_nu_log; % LIKELIHOOD likelihood_log = zeros(1,ndraw); for i = 1:ndraw likelihood_log(i) = ARCHnct_likelihood(Rt,zt,draws(i,:)); end GD_log = source_pdf_log - prior_pdf_log - likelihood_log; constant = max(GD_log); ML_GD = -constant - log(mean(exp(GD_log - constant)))
3.2 Kode bisection function mv = bisection_nu(hpv,zt) % Tujuan : Mencari akar dari F’(v)= 0 menggunakan metode bagi dua % % ------------------------------------------------------------------------% Masukan : hpv = nilai prior v % zt = [z_1, z_2, z_3, ….,z_T] % ---------------------------------------------------------------------
36
% keluaran : mv = akar penyelesaian eps_step = 1e-2; bb = 2.1; ba = 100; if diffFnu(hpv,zt,bb) == 0 % derivatif pertama mv = bb; return; elseif diffFnu(hpv,zt,ba) == 0 mv = ba; return; elseif diffFnu(hpv,zt,bb)*diffFnu(hpv,zt,ba) > 0 error( 'diffFnu(a) and diffFnu(b) do not have opposite signs' ); end while abs(ba - bb) >= eps_step x = (ba + bb)/2; if diffFnu(hpv,zt,x) == 0 mv = x; return; elseif diffFnu(hpv,zt,x)*diffFnu(hpv,zt,ba) < 0 bb = x; else ba = x; end end mv=x;
3.3 Kode bisection NCT function mab = bisection_nct(Rt,R2,a,b,c,hp1,zt,par) % Tujuan : Mencari akar dari F’(a)=0 atau F’(b)=0 menggunakan metode % bagi dua % % --------------------------------------------------------------------% Masukan : Rt = returns % R2 = returns^2 % a = nilai a % b = nilai b % c = nilai miu % hp1 = nilai prior % zt = [z_1, z_2, z_3, ….,z_T] % par = 'a' atau 'b' % ---------------------------------------------------------------------% keluaran : mab = akar penyelesaian eps_step = 1e-2; if par == 'a' bb = 1e-5; ba = 1; elseif par == 'b' bb = 0; ba = 1; end if diffARCH_nct(Rt,R2,a,b,c,hp1,zt,bb,par) == 0 %derivatif pertama mab = bb; return; elseif diffARCH_nct(Rt,R2,a,b,c,hp1,zt,ba,par) == 0 mab = ba; return; elseif diffARCH_nct(Rt,R2,a,b,c,hp1,zt,ba,par)*diffARCH_nct(Rt,R2,a,b,c,hp1,zt,bb,par) > 0 error( 'diffARCH(ba) and diffARCH(bb) do not have opposite signs' );
37
end while abs(bb - ba) >= eps_step x = (ba + bb)/2; if diffARCH_nct(Rt,R2,a,b,c,hp1,zt,c,par) == 0 mab = x; return; elseif diffARCH_nct(Rt,R2,a,b,c,hp1,zt,c,par)*diffARCH_nct(Rt,R2,a,b,c,hp1,zt,ba,par) < 0 bb = x; else ba = x; end end mab = x;
3.4 Kode untuk turunan pertama F() function y = diffFnu(hpv,zt,nu) % Tujuan : Mengitung F’(v) % ----------------------------------------------------------------------% Masukan : hpv = nilai prior % zt = [z_1, z_2, z_3, ….,z_T] % nu = nilai v % ----------------------------------------------------------------------% keluaran : y = nilai turunan pertama T=length(zt); %derivatif pertama y = T/2*(log(nu/2)+1-psi(nu/2))-0.5*sum(log(zt)+1./zt)+(hpv(1)-1)/nu-hpv(2);
4.5 Kode untuk turunan pertama F(a) dan F(b) function Fab = diffARCH_nct(Rt,R2,a,b,c,hp1,zt,bts,par) % Tujuan : Mengitung F’(a) atau F’(b) % --------------------------------------------------------------------% Masukan : Rt = returns % R2 = returns^2 % a = nilai a % c = nilai miu % hp1 = nilai prior % zt = [z_1, z_2, z_3, ….,z_T] % bts = batas kiri/kanan interval pada metode bagi dua % par = 'a' atau 'b' % --------------------------------------------------------------------% keluaran : Fab = nilai turunan pertama if par =='a'% derivatif pertama terhadap a a = bts; lamd = hp1; Fab = -1/(2*a)-0.5*sum(1./(a+b*R2(1:end-1)))... +0.5*(1-b)*R2(1)/(a^2*zt(1))... -2*(1-b)^2*R2(1)./(a^3+zt(1).^2)... +0.5*sum(((R2(2:end))./((a+b*R2(1:end-1)).^2.*zt(2:end)))... -c*Rt(2:end)./(a+b*R2(1:end-1)).^(1.5).*sqrt(zt(2:end)))... -lamd; elseif par == 'b' % derivatif pertama terhadap b b = bts; alp = hp1(1); bet=hp1(2); Fab = -1/(2*(1-b))... -0.5*sum(R2(1:end-1)./(a+b*R2(1:end-1)))... +R2(1)/(2*a*zt(1))... -0.5*c*Rt(1)/sqrt(1-b/a*zt(1))*a*zt(1)... +0.5*sum(R2(1:end-1).*R2(2:end)./((a+b*R2(1:end... -1)).^2.*zt(2:end)))-1.5*sum(R2(1:end... -1).^2*c.*Rt(2:end)./(a+b*R2(1:end-1)).^2.5.*zt(2:end))...
38
+(alp-1)/b-(bet-1)/(1-b); end
3.6 Kode likelihood NCT function [log_likelihood] = ARCHnct_likelihood(R,zt,draws) % Tujuan : Mencari log likelihood ARCH(1) NCT % % --------------------------------------------------------------------% Masukan : R = returns % zt = [z_1, z_2, z_3, ….,z_T] % draws = [av bv cv nuv] % ---------------------------------------------------------------------% keluaran : log_likelihood T = max(size(R)); sig2 = zeros(T,1); sig2(1) = draws(1)/(1-draws(2)); var(1) = sig2(1)*zt(1); y2(1) = (R(1)-sqrt(var(1))*draws(3))^2; log_likelihood = -0.5*log(2*pi*var(1)) -0.5*y2(1)/var(1); for i=2:T sig2(i) = draws(1) + draws(2)*R(i-1)^2; var(i) = sig2(i)*zt(i); y2(i) = (R(i)-sqrt(var(i))*draws(3))^2; log_likelihood = log_likelihood -0.5*log(2*pi*var(i)) -0.5*y2(i)/var(i); end
3.7 Kode untuk memanggil hasil estimasi clc clear all %data data1=csvread('jpy_kursbeli.csv'); data1 = data1(end-1471:end); data2=csvread('eur_kursbeli.csv'); data3 = data3(end-1471:end); Rjpy=100*price2ret(data1); Rjpy=Rjpy-mean(Rjpy); Reur=100*price2ret(data2); Reur=Reur-mean(Reur); lamd=1; alpb=2.5; betb=3; alpnu=16; betnu=0.8; muc=0; varc=1; HP=[lamd alpb betb muc varc alpnu betnu]; hasil vol = zv = av = bv = cv = nuv = draws
= archnct_mcmc(Rjpy,HP); hasil.vol ; hasil.zv; hasil.av; hasil.bv; hasil.cv; hasil.nuv; = hasil.draws;
% ========== Save data ==================================================== dlmwrite('F:\nct\HASIL\eur_nct_draws.csv',draws,'delimiter',',','precision','%0 .16f')
39
3.8 Kode Pendukung Kode truncnormrnd, IACT, momentg, HPD, NSE, mprint, dan mvn_pdf_log dapat dilihat dalam Nugroho (2014).
40
Lampiran 4. Kode Matlab untuk Estimasi Model ARCH(1) Skewed Student-t 4.1 Kode Utama function hasil=archskt_mcmc(Rt,HP) % % % % % % % % % % % %
Tujuan: Mengestimasi parameter-parameter dalam model ARCH(1) R_t = sigma_t*(k*(zt-miu_z)+zt^0.5*xhi_t), xhi_t~N(0,1) sigma_t^2 = a + b*R_{t-1}^2 ------------------------------------------------------------------------Algoritma: MCMC ------------------------------------------------------------------------Masukan: Rt = Returns = 100*[(log(S_t)-log(S_{t-1}))- 1/T*sum((log(S_t)log(S_{t-1}))], dimana S_t = nilai kurs pada saat t HP = Hyperparameter = [lambda alpha_b beta_b miu_k var_k alpha_nu beta_nu], untuk prior a ~ exp(lambda), b ~ beta(alpha_b,beta_b), k ~ Normal(miu_k,var_k), dan nu ~ Gamma(alpha_nu, beta_nu)
% ------------------------------------------------------------------------% Keluaran: hasil.draws = draws % ----- I: Inisialisasi T = length(Rt); R2 = Rt.^2; % prior untuk a lamd = HP(1); % prior untuk b alp = HP(2); bet = HP(3); %prior untuk k muk=HP(4); vark=HP(5); % prior untuk nu alpnu = HP(6); betnu = HP(7); % nilai awal a = 0.1; b = 0.1; nu = 10; zt = 1./gamrnd(nu/2,2/nu,T,1); mz=nu/(nu-2); sig2t=[a/(1-b); a+b*R2(1:T-1)]; % banyak replikasi Nits = 15000; BI = 5000; N = Nits-BI; % alokasi penyimpanan nilai-nilai acak av = zeros(N,1); bv = zeros(N,1); kv = zeros(N,1); nuv = zeros(N,1); zv = zeros(T,1); vol = zeros(T,1); % ----- Algoritma MCMC. Step 1: Membangkitkan Rantai Markov tic for its = 1:Nits % --- pembangkitan sampel k At = sig2t.^(0.5).*(zt-mz);
41
Bt sum1 sum2 Vk = Mk = k =
= sig2t.*zt; = sum(At.^2./Bt); = sum(-1./2*vark); 1./(sum1+sum2); Vk*sum(At.*Rt./Bt+muk/vark); normrnd(Mk,sqrt(Vk),1,1);
% --- pembangkitan sampel z alpz = (nu+1)/2; betz = nu/2+R2/2+k*Rt.*sig2t.^(-0.5)*mz... +k^2.*mz^2/2; % pembangkitan proposal zprop = 1./gamrnd(alpz,1./betz); % mengevaluasi probabilitas penerimaan gzold = k*Rt.*sig2t.^(-0.5)-0.5*k^2*zt+k^2*mz*zt; gznew = k*Rt.*sig2t.^(-0.5)-0.5*k^2*zprop+k^2*mz*zprop; ratio = exp(gznew-gzold); % pembangkitan variabel acak seragam u = rand(T,1); % pembaharuan minrat = min(1,ratio); ind = find(u<=minrat); zt(ind) = zprop(ind); % mencari rata-rata dan variansi untuk proposal bersyarat hpv = [alpnu betnu]; mv = bisection_nuSkt(Rt,sig2t,k,hpv,zt); % rata-rata ddv =sum(-4*k./(sig2t.^(0.5).*zt)*(mv-2)^(-4).*((mv-2)*(Rtsig2t.^(0.5).*k.*zt... +sqrt(sig2t).*k.*(mv/(mv-2))).*sqrt(sig2t).*k) ); d2v = ddv +0.5*T/mv-T/4*psi(1,mv/2)-(alpnu-1)/mv^2; Dv = min(-0.0001,d2v); s2v = -1/Dv; %variansi % algoritma IC-MH % pembangkitan proposal nu* pr = truncnormrnd(1,mv,sqrt(s2v),2.1,40); % mengevaluasi probabilitas penerimaan log_pv = -0.5*sum((2*Rt*k.*sig2t.^(-0.5)*pr/(pr-2)... +k^2.*(-2.*zt*pr/(pr-2)+(pr/(pr-2))^2))./zt)... +0.5*pr*T*log(pr/2)-T*gammaln(pr/2) ... -pr/2*sum(log(zt)+1./zt)+(alpnu-1)*log(pr)-betnu*pr; post_pr = exp(log_pv); log_pv = -0.5*sum((2*Rt*k.*sig2t.^(-0.5)*(nu./(nu-2))... +k^2*(-2*zt*(nu/(nu-2))+(nu./(nu-2))^(2)))./zt)... +0.5*nu*T*log(nu/2)-T*gammaln(nu/2)-nu/2 ... *sum(log(zt)+1./zt)+(alpnu-1)*log(nu)-betnu*nu; post_o = exp(log_pv); ratio = post_pr/post_o; ap = min(1,ratio); % pembangkitan variabel acak seragam u = rand(1); % pembaruan if u <= ap, nu = pr; end mz = nu/(nu-2); % --- pembangkitan sampel a % mencari rata-rata dan variansi untuk distribusi proposal hp1 = lamd; m_a = bisection_skt(Rt,R2,a,b,mz,k,hp1,zt,'a'); % rata-rata d2a = 1/(2*m_a^2)+0.5*sum(1./(m_a+b*R2(1:T-1)).^2)... -0.5*(1-b)*R2(1)/(m_a^3*zt(1))-(0.25*k*Rt(1)... *(1-b)^2)./(((1-b)/m_a*zt(1))^(1.5)*m_a^(4))... +k*Rt(1)*(1-b)./sqrt((1-b)./m_a*m_a^(3)... -(0.75)*k*Rt(1)*mz)./((m_a/(1-b))^(2.5)*zt(1)*(1-b)^(2))...
42
-sum(R2(2:T)./((m_a+b*R2(1:T-1)).^3.*zt(2:T)))... +0.75*sum(k*Rt(2:T)./(m_a+b*Rt(1:T-1).^(2.5)))... -0.75*sum((k*Rt(2:T).*mz)./(m_a... +b*Rt(1:T-1).^(2.5).*zt(2:T)))); Da = min(-0.0001,d2a); s2_a = -1/Da; % variansi % algoritma IC-MH % pembagkitan proposal a* pr = truncnormrnd(1,m_a,sqrt(s2_a),1e-4,1); % mengevaluasi probabilitas penerimaan log_pa = -0.5*log(pr)-0.5*sum(log(pr+b*R2(1:T-1)))-... 0.5*(1-b)*R2(1)/(pr*zt(1))-... +sqrt((1-b)/pr)*k*Rt(1)... -(k*Rt(1)*mz)./(sqrt(pr/(1-b)).*zt(1))... -0.5*sum(R2(2:T)./((pr+b*R2(1:T-1)).*zt(2:T)) ... +(k*Rt(2:T)./(sqrt(pr+b*R2(1:T-1))))... -(k*Rt(2:T)*mz./(sqrt(pr+b*R2(1:T-1))).*zt(2:T)))... -lamd*pr; % F(pr) post_pr = exp(log_pa); log_pa
post_o
= -0.5*log(a)-0.5*sum(log(a+b*R2(1:T-1)))-... 0.5*(1-b)*R2(1)/(a*zt(1))-... +sqrt((1-b)/a)*k*Rt(1)... -(k*Rt(1)*mz)./(sqrt(a/(1-b)).*zt(1))... -0.5*sum(R2(2:T)./((a+b*R2(1:T-1)).*zt(2:T))... +(k*Rt(2:T)./(sqrt(a+b*R2(1:T-1))))... -(k*Rt(2:T)*mz./(sqrt(a+b*R2(1:T-1))).*zt(2:T)))... -lamd*a; % F(a) = exp(log_pa);
ratio = post_pr/post_o; ap = min(1,ratio); % pembangkitan variabel acak seragam u = rand(1); % pembaruan if u <= ap, a = pr; end % --- pembangkitan sampel b % mencari rata-rata dan variansi untuk distribusi proposal mup = alp; Vp = bet; hp1 = [mup Vp]; m_b = bisection_skt(Rt,R2,a,b,mz,k,hp1,zt,'b'); % rata-rata d2b = -0.5/(1-m_b)^2+0.5*sum((R2(1:T-1).^2)./((a+m_b*R2(1:T... -1)).^2))-0.25*k*Rt(1)/((1-m_b/a).^(1.5)*a^2)... +0.25*k*Rt(1)*mz/((1-m_b/a).^(1.5)*a^2.*zt(1))... -sum(((R2(1:T-1).^2).*(R2(2:T)))./((a+m_b*R2(1:T... -1)).^3.*zt(2:T)))+0.75*sum(R2(1:T... -1).^(2)*k.*Rt(2:T)./(a+m_b*Rt(1:T-1).^(2.5)))... -0.75*sum(R2(1:T-1).^(2)*k.*Rt(2:T)*mz./(a+m_b*Rt(1:T... -1).^(2.5).*zt(2:T)))-(alp-1)/(m_b^2)-(bet-1)/(1-m_b)^2; Db = min(-0.0001,d2b); s2_b = -1/Db; % variansi % algoritma IC-MH % pembangkitan proposal b* pr = truncnormrnd(1,m_b,sqrt(s2_b),0,1); % mengevaluasi probabilitas penerimaan log_pb =0.5.*log(1-pr)-0.5.*sum(log(a+pr.*R2(1:T-1)))... -0.5.*(1-pr).*R2(1)./(a*zt(1))... +sqrt((1-pr)/a).*k.*Rt(1)-k.*Rt(1).*mz./(sqrt(a/(1... -pr)).*zt(1))-0.5.*sum(R2(2:T)./(a+pr.*R2(1:T-... 1).*zt(2:T))+k.*Rt(2:T)./sqrt(a+pr.*R2(1:T-1))... -k.*Rt(2:T).*mz./sqrt(a+pr.*R2(1:T-1)).*zt(2:T))... +(alp-1)*log(pr)+(bet-1)*log(1-pr); % F(pr) post_pr = exp(log_pb);
43
log_pb
post_o
= 0.5.*log(1-b)-0.5.*sum(log(a+b.*R2(1:T-1)))... -0.5.*(1-b).*R2(1)./(a*zt(1))... +sqrt((1-b)/a).*k.*Rt(1)-k.*Rt(1).*mz./(sqrt(a/(1... -b)).*zt(1))-0.5.*sum(R2(2:T)./(a+b.*R2(1:T... -1).*zt(2:T))+k.*Rt(2:T)./sqrt(a+b.*R2(1:T-1))... -k.*Rt(2:T).*mz./sqrt(a+b.*R2(1:T-1)).*zt(2:T))... +(alp-1)*log(b)+(bet-1)*log(1-b); %F(b) = exp(log_pb);
ratio ap
= post_pr/post_o; = min(1,ratio);
% pembangkitan variabel acak seragam u = rand(1); % pembaruan if u <= ap, b = pr; end % pengestimasian sigma volt = [a/(1-b); a+b*R2(1:T-1)]; sig2t = volt; % simpan a, b, k, dan nu if its > BI av(its-BI,1) = a; bv(its-BI,1) = b; kv(its-BI,1) = k; nuv(its-BI,1) = nu; zv = ((its-BI-1)*zv+zt)/(its-BI); vol = ((its-BI-1)*vol+volt)/(its-BI); end end toc % ----draws = MP = SP =
Algoritma MCMC. Step 2: Menghitung rata-rata Monte Carlo [av bv kv nuv]; mean(draws); std(draws);
% ===== Integrated Autocorrelation Time (IACT) ============================ % Berapa banyak sampel yang harus dibangkitkan untuk mendapatkan sampel % yang bebas (seberapa cepat konvergensi simulasi) resultsIAT = IACT(draws); IAT = [resultsIAT.iact]; % ===== Uji Konvergensi Geweke ============================================ idraw1 = round(.1*N); resultCV = momentg(draws(1:idraw1,:)); meansa = [resultCV.pmean]; nsea = [resultCV.nse1]; idraw2 resultCV meansb nseb
= = = =
round(.5*N)+1; momentg(draws(idraw2:N,:)); [resultCV.pmean]; [resultCV.nse1];
CD = (meansa - meansb)./sqrt(nsea+nseb); onetail = 1-normcdf(abs(CD),0,1); pV = 2*onetail; % ===== 95% Highest Posterior Density (HPD) Interval ====================== resultsHPD = HPD(draws,0.05); LB = [resultsHPD.LB]; UB = [resultsHPD.UB]; % % ===== Numerical Standard Error (NSE) ====================================
44
resultsNSE = NSE(draws); NSEd = [resultsNSE.nse]; %====================== Pengaturan Pencetakan Hasil ======================= %----- Statistik Parameter in.cnames = char('a','b','k','nu'); in.rnames = char('Parameter','Mean','SD','LB','UB','IACT','NSE','G-CD','pValue'); in.fmt = '%16.6f'; tmp = [MP; SP; LB; UB; IAT; NSEd; CD; pV]; fprintf(1,'Estimasi menggunakan MCMC dan Uji Diagnostik\n'); mprint(tmp,in); hasil.vol hasil.zv hasil.av hasil.bv hasil.kv hasil.nuv hasil.draws
= = = = = = =
vol; zv; av; bv; kv; nuv; draws;
4.2 Kode bisection function mv = bisection_nuSkt(Rt,sig2t,k,hpv,zt) % Tujuan : Mencari akar dari F’(v)= 0 menggunakan metode bagi dua % % Masukan : Rt = returns % Sig2t = sigma_t^2 % k = nilai parameter skewness k % hpv = nilai prior v % zt = [z_1, z_2, z_3, ….,z_T] % --------------------------------------------------------------------% keluaran : mv = akar penyelesaian eps_step = 1e-2; bb = 2.1; ba = 100; if diffFnuSkt(Rt,sig2t,k,hpv,zt,bb) == 0 % derivatif pertama mv = bb; return; elseif diffFnuSkt(Rt,sig2t,k,hpv,zt,ba) == 0 mv = ba; return; elseif diffFnuSkt(Rt,sig2t,k,hpv,zt,bb).*diffFnuSkt(Rt,sig2t,k,hpv,zt,ba) > 0 error( 'diffFnuSkt(a) and diffFnuSkt(b) do not have opposite signs' ); end while abs(ba - bb) >= eps_step x = (ba + bb)/2; if diffFnuSkt(Rt,sig2t,k,hpv,zt,x) == 0 mv = x; return; elseif diffFnuSkt(Rt,sig2t,k,hpv,zt,x).*diffFnuSkt(Rt,sig2t,k,hpv,zt,ba) < 0 bb = x; else ba = x; end end mv=x;
45
4.3 Kode bisection SKT function mab = bisection_skt(Rt,R2,a,b,mz,k,hp1,zt,par) % Tujuan : Mencari akar dari F’(a)=0 atau F’(b)=0 menggunakan metode % bagi dua % % --------------------------------------------------------------------% Masukan : Rt = returns % R2 = returns^2 % a = nilai a % b = nilai b % mz = nilai miu_z % k = nilai parameter skewness k % hp1 = nilai prior % zt = [z_1, z_2, z_3, ….,z_T] % par = 'a' atau 'b' % ---------------------------------------------------------------------% keluaran : mab = akar penyelesaian eps_step = 1e-2; if par == 'a' bb = 1e-5; ba = 1; elseif par == 'b' bb = 0; ba = 1; end if diffARCH_skt(Rt,R2,a,b,mz,k,hp1,zt,bb,par) == 0 %derivatif pertama mab = bb; return; elseif diffARCH_skt(Rt,R2,a,b,mz,k,hp1,zt,ba,par) == 0 mab = ba; return; elseif diffARCH_skt(Rt,R2,a,b,mz,k,hp1,zt,ba,par)*diffARCH_skt(Rt,R2,a,b,mz,k,hp1,zt,b b,par) > 0 error( 'diffARCH(ba) and diffARCH(bb) do not have opposite signs' ); end while abs(bb - ba) >= eps_step x = (ba + bb)/2; if diffARCH_skt(Rt,R2,a,b,mz,k,hp1,zt,x,par) == 0 mab = x; return; elseif diffARCH_skt(Rt,R2,a,b,mz,k,hp1,zt,x,par)*diffARCH_skt(Rt,R2,a,b,mz,k,hp1,zt,ba ,par) < 0 bb = x; else ba = x; end end mab = x;
46
4.4 Kode untuk turunan pertama F() function y = diffFnuSkt(Rt,sig2t,k,hpv,zt,nu) % Tujuan : Mengitung F’(v) % % ----------------------------------------------------------------------% Masukan : Rt = returns % sig2t = sigma_t^2 % k = nilai parameter skewness k % hpv = nilai prior % zt = [z_1, z_2, z_3, ….,z_T] % nu = nilai v % ----------------------------------------------------------------------% keluaran : y = nilai turunan pertama T = length(zt); y
= sum(2./(sig2t.^(0.5).*zt)*k*(nu-2)^(-2).*(Rt... -sig2t.^(0.5)*k.*zt+sig2t*k*nu/(nu-2)))... +T/2*(log(nu/2)+1-psi(nu/2))-0.5*sum(log(zt)+1./zt)... +(hpv(1)-1)/nu-hpv(2);
4.5 Kode untuk turunan pertama F(a) dan F(b) function Fab = diffARCH_skt(Rt,R2,a,b,mz,k,hp1,zt,bts,par) % Tujuan : Mengitung F’(a) atau F’(b) % --------------------------------------------------------------------% Masukan : Rt = returns % R2 = returns^2 % a = nilai a % b = nilai b % mz = nilai miu_z % k = nilai parameter skewness k % hp1 = nilai prior % zt = [z_1, z_2, z_3, ….,z_T] % bts = batas kiri/kanan interval pada metode bagi dua % par = 'a' atau 'b' % --------------------------------------------------------------------% keluaran : Fab = nilai turunan pertama if par =='a'% derivatif pertama terhadap a a = bts; lamd = hp1; Fab = -1/(2*a)-0.5*sum(1./(a+b*R2(1:end-1)))... +0.5*(1-b)*R2(1)/(a^2*zt(1))... -0.5*k.*(1-b)^2*Rt(1)./((sqrt(1-b)./a).*a.^(2))... +0.5*k.*Rt(1).*mz./sqrt(a./(1-b)).*zt(1).*(1-b)... +0.5*sum(((R2(2:end))./((a+b*R2(1:end-1)).^2.*zt(2:end)))... -0.5*k*Rt(2:end)./(a+b*R2(1:end-1)).^(1.5)... +0.5*k.*Rt(2:end).*mz./(a+b*R2(1:end... -1)).^(1.5).*zt(2:end))-lamd; elseif par == 'b' % derivatif pertama terhadap b b = bts; alp = hp1(1); bet=hp1(2); Fab = -1/(2*(1-b))... -0.5*sum(R2(1:end-1)./(a+b*R2(1:end-1)))... +R2(1)/(2*a*zt(1))... -0.5*k*Rt(1)/sqrt((1-b)/a).*a... +0.5*k*Rt(1).*mz./(sqrt((1-b)./a).*a*zt(1))... +0.5*sum(R2(1:end-1).*R2(2:end)./((a+b*R2(1:end... -1)).^2.*zt(2:end)))... -0.5*sum(R2(1:end-1)*k.*Rt(2:end)./(a+b*R2(1:end-1)).^1.5)... +0.5*sum(R2(1:end-1)*k.*Rt(2:end)./(a+b*R2(1:end... -1)).^1.5.*zt(2:end))+(alp-1)/b-(bet-1)/(1-b);
47
end
4.6 Kode untuk memanggil hasil estimasi clc clear all %data data1=csvread('jpy_kursbeli.csv'); data1 = data1(end-1471:end); data2=csvread('eur_kursbeli.csv'); data3 = data3(end-1471:end); Rjpy=100*price2ret(data1); Rjpy=Rjpy-mean(Rjpy); Reur=100*price2ret(data2); Reur=Reur-mean(Reur); lamd=1; alpb=2.5; betb=3; alpnu=16; betnu=0.8; muk=0; vark=1; HP=[lamd alpb betb muk vark alpnu betnu]; hasil = archskt_mcmc(Reur,HP); draws = hasil.draws; % ========== Save data ==================================================== dlmwrite('D:\SKRIPSI\skt\HASIL\jpy_skt_draws.csv',draws,'delimiter',',','precis ion','%0.16f')
4.7 Kode Pendukung Kode truncnormrnd, IACT, momentg, HPD, NSE, dan mprint dapat dilihat dalam Nugroho (2014).
48
Lampiran 5. Sertifikat Seminar
49
50
