new; library pgraph optmum; pqgwin many; cls; t = 156; load data[t,1] = c:\data\real_int.txt; @ real intrest rates, 1974:1 ~ 1986:12@ rint = data; ymat = rint; t = rows(ymat); s = 2; @ # of regimes @ //xy(t,rint); stop; /*=============== Estimation ===============*/ @ Deriving desirable initial values using OLS @ mu1_in = -3; mu2_in = 4; phi1_in = 0.1; phi2_in = 0.1; sigma1_in = ln(5); sigma2_in = ln(5); p11_in = 3; @ Pr(S(t)=1|S(t-1)=1) @ p22_in = 3; @ Pr(S(t)=2|S(t-1)=2) @ prmtr_in = mu1_in|mu2_in|phi1_in|phi2_in|sigma1_in|sigma2_in|p11_in|p22_in; @ initial values @ {xout,fout,gout,cout}=optmum(&lik,prmtr_in); xout_fnl = trans(xout); cov = inv(hessp(&lik,xout)); grad = gradfd(&trans,xout); cov_fnl = grad*cov*grad'; stde = sqrt(diag(cov_fnl)); t_val = xout_fnl./stde; /*===============Filtering ===================*/ Prbmat = filter(xout); @ Filtered probability @ seq = seqa(74+1/12,1/12,t); graphset; _ptitlht = 0.25; begwind; window(2,1,0); setwind(1); title(" St = 1 "); xy(seq,prbmat[.,1]); setwind(2); title(" St = 2 "); xy(seq,prbmat[.,2]); endwind; /*==========Output===========*/ "================================================"; " Estimates S.D. t- values"; "------------------------------------------------"; xout_fnl~stde~t_val ; "------------------------------------------------"; end; /*========== Likelihood function ======*/ proc lik(prmtr0); local prmtr,mu1,phi1,mu2,phi2,sigma1,sigma2, p11,p22,p21,p12, PPR_TR,AA,PPROB_T,EN,prb_LL1,prb_LL2, lnL,itr,prb_tL11,prb_tL12,prb_tL21,prb_tL22, et11,et12,et21,et22,pdf11,pdf21,pdf12,pdf22, jpdf11,jpdf21,jpdf12,jpdf22,pdf, lnf,prb_tt1,prb_tt2; prmtr = trans(prmtr0); mu1 = prmtr[1]; @ unconditional mean in state 1 @ mu2 = prmtr[2]; @ unconditional mean in state 2 @ phi1 = prmtr[3]; @ AR coefficient in state 1 @ phi2 = prmtr[4]; @ AR coefficient in state 2 @ sigma1 = prmtr[5]; @ variance in state 1 @ sigma2 = prmtr[6]; @ variance in state 2 @ p11 = prmtr[7]; @ Pr(S(t)=1|S(t-1)=1) @ p22 = prmtr[8]; @ Pr(S(t)=2|S(t-1)=2) @ p12 = 1 - p11; p21 = 1 - p22; PPR_TR = p11~p12| p21~p22; @ s by s @ AA = (eye(s)-PPR_TR)|ones(1,s); @ (s+1) by s @ EN = zeros(s,1)|1; PPROB_T = inv(AA'AA)*AA'en; @ Unconditional(initial) Probability, s by 1 @ prb_LL1 = PPROB_T[1]; @ Pr[S(t-1)=1|I(t-1)] @ prb_LL2 = PPROB_T[2]; @ Pr[S(t-1)=2|I(t-1)] @ lnL = 0; itr = 2; do until itr>t; /* Pr[S(t),S(t-1)|I(t-1)] */ prb_tL11 = prb_LL1*p11;@ Pr[S(t)|S(t-1)] times Pr[S(t-1)| I(t-1) ] @ prb_tL12 = prb_LL1*p12; prb_tL21 = prb_LL2*p21; prb_tL22 = prb_LL2*p22; /* prediction error conditioned on S(t) and S(t-1) */ et11 = ymat[itr] - mu1 - phi1*(ymat[itr-1] - mu1); et12 = ymat[itr] - mu2 - phi2*(ymat[itr-1] - mu1); et21 = ymat[itr] - mu1 - phi1*(ymat[itr-1] - mu2); et22 = ymat[itr] - mu2 - phi2*(ymat[itr-1] - mu2); /* f[y|S(t),S(t-1),I(t-1)]: conditional density */ pdf11 = exp(lnpdfmvn(et11',sigma1)); pdf12 = exp(lnpdfmvn(et12',sigma2)); pdf21 = exp(lnpdfmvn(et21',sigma1)); pdf22 = exp(lnpdfmvn(et22',sigma2)); /* Pr[y(t),S(t),S(t-1)|I(t-1)]: joint conditional density */ jpdf11 = pdf11*prb_tL11; jpdf12 = pdf12*prb_tL12; jpdf21 = pdf21*prb_tL21; jpdf22 = pdf22*prb_tL22; /* Probability Density Function for y(t) */ pdf = jpdf11 + jpdf12 + jpdf21 + jpdf22; lnf = ln(pdf); @ log density @ lnL = lnL + lnf; /* Pr[S(t)|I(t)] */ prb_tt1 = (jpdf11 + jpdf21)/pdf; prb_tt2 = (jpdf12 + jpdf22)/pdf; prb_LL1 = prb_tt1; prb_LL2 = prb_tt2; itr = itr + 1; endo; retp(-lnL); endp; /*========== Procedure 1:Defining Likelihood function ======*/ proc Filter(prmtr0); local prmtr,mu1,phi1,mu2,phi2,sigma1,sigma2, p11,p22,p21,p12, PPR_TR,AA,PPROB_T,EN,prb_LL1,prb_LL2, lnL,itr,prb_tL11,prb_tL12,prb_tL21,prb_tL22, et11,et12,et21,et22,pdf11,pdf21,pdf12,pdf22,pdf, jpdf11,jpdf21,jpdf12,jpdf22, lnf,prb_tt1,prb_tt2,Prob_tt; prmtr = trans(prmtr0); mu1 = prmtr[1]; @ unconditional mean in state 1 @ mu2 = prmtr[2]; @ unconditional mean in state 2 @ phi1 = prmtr[3]; @ AR coefficient in state 1 @ phi2 = prmtr[4]; @ AR coefficient in state 2 @ sigma1 = prmtr[5]; @ variance in state 1 @ sigma2 = prmtr[6]; @ variance in state 2 @ p11 = prmtr[7]; @ Pr(S(t)=1|S(t-1)=1) @ p22 = prmtr[8]; @ Pr(S(t)=2|S(t-1)=2) @ p12 = 1 - p11; p21 = 1 - p22; PPR_TR = p11~p12| p21~p22; @ s by s @ AA = (eye(s)-PPR_TR)|ones(1,s); @ (s+1) by s @ EN = zeros(s,1)|1; PPROB_T = inv(AA'AA)*AA'en; @ Unconditional Probability, s by 1 @ prb_LL1 = PPROB_T[1]; @ Pr[S(t-1)=1|I(t-1)] @ prb_LL2 = PPROB_T[2]; @ Pr[S(t-1)=2|I(t-1)] @ Prob_tt = zeros(t,s); @ storage for filtered prb @ Prob_tt[1,1] = prb_LL1; Prob_tt[1,2] = prb_LL2; itr = 2; do until itr>t; /* Pr[S(t),S(t-1)|I(t-1)] */ prb_tL11 = prb_LL1*p11;@ Pr[S(t)|S(t-1)] times Pr[S(t-1)| I(t-1) ] @ prb_tL12 = prb_LL1*p12; prb_tL21 = prb_LL2*p21; prb_tL22 = prb_LL2*p22; /* prediction error conditioned on S(t) and S(t-1) */ et11 = ymat[itr] - mu1 - phi1*(ymat[itr-1] - mu1); et12 = ymat[itr] - mu2 - phi2*(ymat[itr-1] - mu1); et21 = ymat[itr] - mu1 - phi1*(ymat[itr-1] - mu2); et22 = ymat[itr] - mu2 - phi2*(ymat[itr-1] - mu2); pdf11 = exp(lnpdfmvn(et11',sigma1)); pdf12 = exp(lnpdfmvn(et12',sigma2)); pdf21 = exp(lnpdfmvn(et21',sigma1)); pdf22 = exp(lnpdfmvn(et22',sigma2)); jpdf11 = pdf11*prb_tL11; jpdf12 = pdf12*prb_tL12; jpdf21 = pdf21*prb_tL21; jpdf22 = pdf22*prb_tL22; /* Probability Density Function for ymat[itr] */ pdf = jpdf11 + jpdf12 + jpdf21 + jpdf22; /* Pr[S(t)|I(t)] */ prb_tt1 = (jpdf11 + jpdf21)/pdf; prb_tt2 = (jpdf12 + jpdf22)/pdf; Prob_tt[itr,1] = prb_tt1; Prob_tt[itr,2] = prb_tt2; prb_LL1 = prb_tt1; prb_LL2 = prb_tt2; itr = itr + 1; endo; retp(Prob_tt); endp; /*======== Trans ========*/ proc trans(prmtr1); local prmtr; prmtr = prmtr1; prmtr[3:4] = prmtr1[3:4]./ (abs(prmtr1[3:4])+ 1); @ phi @ prmtr[5:6] = exp(prmtr1[5:6]); @ sigma @ prmtr[7:8] = exp(prmtr1[7:8])./ (exp(prmtr1[7:8])+ 1); @ prob @ retp(prmtr); endp;