/*========= Changes in U.S. Inflation Persistence =========== This program is designed for 2 permanant structural changes model. Written by Kyu Ho Kang, March 2008 data : US GDP deflator inflation < Model > Yt= Z_t + X_t Z_t = mu_S_t + Z_t-1 + Vt, Vt ~ N(0, sig2v_S(t)) X_t = b1_St*X_t-1+ b2_St*X_t-2 + Wt, Wt ~ N(0,sig2w_St) Vt and Wt are correlated. St = { 1 ,2 ,3} Prb[ S(t)=1 | S(t-1)=1] = p11 Prb[ S(t)=2 | S(t-1)=2] = p22 ==================================================*/ new; cls; library pgraph optmum; pqgwin many; start=1; cut=47; k=3;@ dimension of Beta @ n=1;@ dimension of Y_t @ pv=1.96; j=40; t=237; load data3[t,1] = us_def_inf.txt; ymat=(data3[cut+1:t]); t=rows(ymat); xy(t,ymat); prmtr_in={ 0.20151076 0.075978685 -0.11315029 -0.26992024 -0.56722643 -0.30285285 -2.9880988 0.27118817 -2.4716655 -1.3472776 -1.7647866 -1.5925288 3.8238695 4.0208352 3.949452 3.6968571 3.7301786 0.099672134 -0.040651884 0.0016778200 }; {xout, fout,gout, cout }=optmum(&lik_f, prmtr_in); xout_fnl=trans(xout); cov=inv(hessp(&lik_f,xout)); grad=gradp(&trans,xout); cov_fnl=grad*cov*grad'; stde=sqrt(diag(cov_fnl)); t_val=xout_fnl./(stde); p11=xout_fnl[13]; p22=xout_fnl[14]; duration1=duration(p11); grad_d=gradp(&duration,p11); cov_d=grad_d*cov_fnl[13,13]*grad_d'; stde_d1=sqrt(diag(cov_d)); duration2=duration(p22); grad_d=gradp(&duration,p22); cov_d=grad_d*cov_fnl[14,14]*grad_d'; stde_d2=sqrt(diag(cov_d)); prmtr=xout_fnl; aa1=prmtr[1]; aa2=prmtr[2]; aa3=prmtr[3]; bb1=prmtr[4]; bb2=prmtr[5]; bb3=prmtr[6]; FF1=(aa1~bb1)| ( 1 ~ 0 ); Psi1=inv(eye(k-1)-FF1); Psi1=Psi1[1,1]; FF2=(aa2~bb2)| ( 1 ~ 0 ); Psi2=inv(eye(k-1)-FF2); Psi2=Psi2[1,1]; FF3=(aa3~bb3)| ( 1 ~ 0 ); Psi3=inv(eye(k-1)-FF3); Psi3=Psi3[1,1]; /*==== Test for variance in transitory component ============*/ prmtr_in0={ 0.20479237 0.065543219 -0.10818384 -0.28234788 -0.52564530 -0.29044911 -2.8110621 0.19169496 -2.5229477 -1.5275590 3.8140202 4.0284334 5.9463767 5.5300149 5.8599186 0.10089309 -0.039991183 0.0014164169 }; {xout0, fout0,gout0, cout0 }=optmum(&lik_f0, prmtr_in0); LR_stats=2*(fout0-fout); prmtr_in1={ 0.45913752 -0.65464983 -0.14221878 0.078882753 -0.60326841 -0.27492810 -2.2732085 -0.089447637 -190.13671 -1.5757578 4.4981045 1.8868541 1601.7306 50.529080 77.951312 0.10084798 -0.68305903 -5.8936623e-005 }; {xout1, fout1,gout1, cout1 }=optmum(&lik_f1, prmtr_in1); LR_stats1=2*(fout1-fout); cls; "======================================"; " Estimated values S.D."; "======================================"; xout_fnl~stde~t_val; "======================================"; "lnL";;-fout; "======================================"; "Expected Duration of Regime 1";;duration1; "Expected Duration of Regime 2";;duration2; "======================================"; "The 1st Breakpoint ";;47.25+(cut+duration1)*0.25; "The 2st Breakpoint ";;47.25+(cut+duration1+duration2)*0.25; "======================================"; "Ratio of Standard deviation(Permanent/Temporary)"; "Regime 1";; xout_fnl[7]/xout_fnl[10]; "Regime 2";; xout_fnl[8]/xout_fnl[11]; "Regime 3";; xout_fnl[9]/xout_fnl[12]; "======================================"; "Sum of AR coeficients"; "Regime 1";;xout_fnl[1]+xout_fnl[4]; "Regime 2";;xout_fnl[2]+xout_fnl[5]; "Regime 3";;xout_fnl[3]+xout_fnl[6]; "======================================"; "Long Run Response of Temporary shock"; "Regime 1";;Psi1; "Regime 2";;Psi2; "Regime 3";;Psi3; "======================================"; "Test H0:sig2w1=sig2w2=sig2w3"; "LR statistics";;abs(LR_stats); "p value";;cdfchic(abs(LR_stats),2); "======================================"; "Test H0:sig2v1=sig2v2=sig2v3"; "LR statistics";;abs(LR_stats1); "p value";;cdfchic(abs(LR_stats1),2); "======================================"; "BIC";;fout+0.5*rows(xout_fnl)*ln(T); "AIC";;fout+rows(xout_fnl); "======================================"; @========== Model Transformation into ARIMA =============@ arma_fnl=arma(xout_fnl); grad_am=gradp(&arma,xout_fnl); cov_am=grad_am*cov_fnl*grad_am'; stde_am=sqrt(diag(cov_am)); t_val_am=arma_fnl./stde_am; @ St=1 @ phi11=arma_fnl[1]; phi21=arma_fnl[2]; theta11=arma_fnl[3]; theta21=arma_fnl[4]; sig2e1=arma_fnl[5]; F1=(phi11~phi21~theta11~theta21)| (1~0~0~0)| zeros(1,4)| (0~0~1~0); H1=1|0|1|0; Psi1=inv(eye(4)-F1)*H1; Psi1=Psi1[1,1]; @ Long-Run effect @ @ St=2 @ phi12=arma_fnl[6]; phi22=arma_fnl[7]; theta12=arma_fnl[8]; theta22=arma_fnl[9]; sig2e2=arma_fnl[10]; F2=(phi12~phi22~theta12~theta22)| (1~0~0~0)| zeros(1,4)| (0~0~1~0); H2=1|0|1|0; imp=0; tmpt=H2; i=1; do until i>j; imp=imp+tmpt[1,1]; tmpt=F2*tmpt; i=i+1; endo; Psi2=inv(eye(4)-F2)*H2; Psi2=Psi2[1,1]; @ Long-Run effect @ @ St=3 @ phi13=arma_fnl[11]; phi23=arma_fnl[12]; theta13=arma_fnl[13]; theta23=arma_fnl[14]; sig2e3=arma_fnl[15]; F3=(phi13~phi23~theta13~theta23)| (1~0~0~0)| zeros(1,4)| (0~0~1~0); H3=1|0|1|0; Psi3=inv(eye(4)-F3)*H3; Psi3=Psi3[1,1]; @ Long-Run effect @ @ ================= Impulse Response ================== @ @ Regime 1 @ itr=1; tmp=H1; cml_imp1=zeros(j,3); imp1=zeros(j,1); do until itr>j; imp=square1(xout); @====================@ proc square1(input); local out,FF,i,FFF,H,tmpt,imp,am_input,phi_hat; input=trans(input); am_input=arma(input); phi_hat=am_input[1:4]; FF=phi_hat'| (1~0~0~0)| zeros(1,4)| (0~0~1~0); H=1|0|1|0; imp=0; tmpt=H; i=1; do until i>itr; imp=imp+tmpt[1,1]; tmpt=FF*tmpt; i=i+1; endo; retp(imp); endp; @=====================@ cml_imp1[itr,2]=imp; grad=gradp(&square1,xout); cov_imp=grad*cov*grad'; @ Delta method @ stde_imp=sqrt(diag(cov_imp)); cml_imp1[itr,1]=imp-pv*stde_imp[1]; cml_imp1[itr,3]=imp+pv*stde_imp[1]; itr=itr+1; endo; @ Regime 2 @ itr=1; cml_imp2=zeros(j,3); do until itr>j; imp=square2(xout); @====================@ proc square2(input); local out,FF,i,FFF,H,tmpt,imp,am_input,phi_hat; input=trans(input); am_input=arma(input); phi_hat=am_input[6:9]; FF=phi_hat'| (1~0~0~0)| zeros(1,4)| (0~0~1~0); H=1|0|1|0; imp=0; tmpt=H; i=1; do until i>itr; imp=imp+tmpt[1,1]; tmpt=FF*tmpt; i=i+1; endo; retp(imp); endp; @=====================@ cml_imp2[itr,2]=imp; grad=gradp(&square2,xout); cov_imp=grad*cov*grad'; @ Delta method @ stde_imp=sqrt(diag(cov_imp)); cml_imp2[itr,1]=imp-pv*stde_imp[1]; cml_imp2[itr,3]=imp+pv*stde_imp[1]; itr=itr+1; endo; @ Regime 3 @ itr=1; cml_imp3=zeros(j,3); do until itr>j; imp=square3(xout); @====================@ proc square3(input); local out,FF,i,FFF,H,tmpt,imp,am_input,phi_hat; input=trans(input); am_input=arma(input); phi_hat=am_input[11:14]; FF=phi_hat'| (1~0~0~0)| zeros(1,4)| (0~0~1~0); H=1|0|1|0; imp=0; tmpt=H; i=1; do until i>itr; imp=imp+tmpt[1,1]; tmpt=FF*tmpt; i=i+1; endo; retp(imp); endp; @=====================@ cml_imp3[itr,2]=imp; grad=gradp(&square3,xout); cov_imp=grad*cov*grad'; @ Delta method @ stde_imp=sqrt(diag(cov_imp)); cml_imp3[itr,1]=imp-pv*stde_imp[1]; cml_imp3[itr,3]=imp+pv*stde_imp[1]; itr=itr+1; endo; tim=seqa(0,1,j); cml_imp=cml_imp1~cml_imp2~cml_imp3; graphset; _ptitlht=0.25; begwind; window(2,2,0); setwind(1); title(" Impulse 1 "); xy(tim,cml_imp1~zeros(j,1)); setwind(2); title(" Impulse 2 "); xy(tim,cml_imp2~zeros(j,1)); setwind(3); title(" Impulse 3 "); xy(tim,cml_imp3~zeros(j,1)); setwind(4); _plctrl={1,1,1,1}; _pltype={6,6,6,6}; _plegctl={2,0.5,6,3}; _plegstr="Regime 1\0"\ "Regime 2\0"\ "Regime 3\0"\ "Zero line"; title(" Comparison "); xy(tim,cml_imp1[.,2]~cml_imp2[.,2]~cml_imp3[.,2]~zeros(j,1)); endwind; "======================================"; " ARIMA Model"; " Estimated values S.D. t values"; "======================================"; arma_fnl[1:5]~stde_am[1:5]~t_val_am[1:5]; arma_fnl[6:10]~stde_am[6:10]~t_val_am[6:10]; arma_fnl[11:15]~stde_am[11:15]~t_val_am[11:15]; "======================================"; "Long run Multiplier at Regime 1";;Psi1~(cml_imp1[j,3]-cml_imp1[j,2])/1.96; "Long run Multiplier at Regime 2";;Psi2~(cml_imp2[j,3]-cml_imp2[j,2])/1.96; "Long run Multiplier at Regime 3";;Psi3~(cml_imp3[j,3]-cml_imp3[j,2])/1.96; "======================================"; cml_imp1[j,2]~cml_imp2[j,2]~cml_imp3[j,2]; "======================================"; /*========================= Kim's Filtering =============================*/ {trend,cycle,ss2mat,ss3mat,SD_resid,eta_mat,f_mat}=filter(xout); p11=xout_fnl[13]; p22=xout_fnl[14]; p12=1-p11; p23=1-p22; p33=1; /*========================= Kim's Smoothing =============================*/ t=rows(ss2mat); prb_h=(1-ss2mat[.,2]-ss3mat[.,2])~ss2mat[.,2]~ss3mat[.,2]; @ Hamilton filtering @ prb_tl=(1-ss2mat[.,1]-ss3mat[.,1])~ss2mat[.,1]~ss3mat[.,1]; prb_k=KIM(prb_h,prb_tl); prb_k3=prb_k[.,2]; @ Pr[St=3 | I(T) ] @ prb_k2=prb_k[.,1]; @ Pr[St=2 | I(T) ] @ prb_k1=1-prb_k2-prb_k3; @ Pr[St=1 | I(T) ] @ prb_k=prb_k1~prb_k2~prb_k3; @ Probability of S.B. at time t @ prb_sb1=zeros(t,1); prb_sb2=zeros(t,1); conf1=zeros(2,1); conf2=zeros(2,1); itr=2; do until itr>t; prb_sb1[itr]=maxc((prb_k2[itr]-prb_k2[itr-1])|0); if sumc(prb_sb1[1:itr-1])<0.025 and sumc(prb_sb1[1:itr])>=0.025; conf1[1]=1947.25+(cut+itr-1)*0.25; endif; if sumc(prb_sb1[1:itr-1])<0.975 and sumc(prb_sb1[1:itr])>=0.975; conf1[2]=1947.25+(cut+itr)*0.25; endif; prb_sb2[itr]=maxc((prb_k3[itr]-prb_k3[itr-1])|0); if sumc(prb_sb2[1:itr-1])<0.025 and sumc(prb_sb2[1:itr])>=0.025; conf2[1]=1947.25+(cut+itr-1)*0.25; endif; if sumc(prb_sb2[1:itr-1])<0.975 and sumc(prb_sb2[1:itr])>=0.975; conf2[2]=1947.25+(cut+itr)*0.25; endif; itr=itr+1; endo; "-----------------------------------------------------------------------------------------------------------------------------"; "95% confidence band for 1st structural break date";;conf1'; "95% confidence band for 2nd structural break date";;conf2'; "-----------------------------------------------------------------------------------------------------------------------------"; prb_sb=prb_sb1~prb_sb2; /*=============== Diagnostic Check =============*/ t0=rows(ss2mat); lag=20; R_sq=SD_resid.*SD_resid;@ Square of Residuals @ gam_mat=zeros(lag,1); gam2_mat=zeros(lag,1); q_mat=zeros(lag,1); @ Ljung-Box Q-statisic for autocorrelation @ p_mat=zeros(lag,1); @ p value @ q2_mat=zeros(lag,1); @ Ljung-Box Q-statisics for ARCH @ p2_mat=zeros(lag,1); @ p value @ jj=1; do until jj>lag; e_j=trim(SD_resid,0,jj) -meanc(trim(SD_resid,0,jj)); e_0=trim(SD_resid,jj,0)-meanc(trim(SD_resid,jj,0)); e2_j=trim(R_sq,0,jj)-meanc(trim(R_sq,0,jj)) ; e2_0=trim(R_sq,jj,0)-meanc(trim(R_sq,jj,0)); gam_j=inv(e_j'e_j)*e_j'e_0; gam2_j=inv(e2_j'e2_j)*e2_j'e2_0; gam_mat[jj]=(gam_j^2)/(t0@-jj@); gam2_mat[jj]=(gam2_j^2)/(t0@-jj@); q_mat[jj]=t0*(t0+2)*sumc(gam_mat[1:jj]); @ Q-statistics @ p_mat[jj]=cdfchic(q_mat[jj],jj); q2_mat[jj]=t0*(t0+2)*sumc(gam2_mat[1:jj]);@ Q-statistics @ p2_mat[jj]=cdfchic(q2_mat[jj],jj); jj=jj+1; endo; "-------------------------------------------------------------------------------------------------------"; " Q for AC p value Q for ARCH p value"; "------------------------------------------------------------------------------------------------------"; q_mat~p_mat~q2_mat~p2_mat; "==================================================="; /*======= Plotting ========*/ t0=rows(ss2mat); tim=seqa(47.25+(cut)*0.25,0.25,t0); ymat=ymat[1:t0]; graphset; _ptitlht=0.25; begwind; window(2,2,0); setwind(1); title(" Prb on St =2"); xy(tim , prb_k2~ss2mat[.,2]~ ones(t0,1)~zeros(t0,1)); setwind(2); title("Stochastic Trend"); xy(tim,trend~ymat); setwind(3); title(" Prb on St=2 and Y"); xy(tim,(prb_k2)*10~ymat); setwind(4); title("Temporary Term"); xy(tim,cycle~zeros(t0,1)); endwind; graphset; _ptitlht=0.25; begwind; window(2,2,0); setwind(1); title(" Prediction error"); xy(tim , eta_mat~zeros(t0,1)); setwind(2); title(" S.D of prediction error"); xy(tim,sqrt(f_mat)); setwind(3); title("Standardized Residuals"); xy(tim,SD_resid~zeros(t0,1)); setwind(4); endwind; graphset; _ptitlht=0.25; begwind; window(2,2,0); setwind(1); title(" Prb on St =1"); xy(tim , (1-prb_k2-prb_k3)~(1-ss2mat[.,2]-ss3mat[.,2])); setwind(2); title(" Prb on St =2"); xy(tim , prb_k2~ss2mat[.,2]); setwind(3); title(" Prb on St=3 "); xy(tim , prb_k3~ss3mat[.,2]); setwind(4); title("Prb on S.B."); xy(tim,prb_sb1*100~prb_sb2*100~ymat); setwind(4); endwind; end; /*================= Lik_f ===================*/ proc lik_f(prmtr1); local prmtr,mmu0,mmu1,aa3,aa1,aa2,ssig2v2,ssig2v3,ssig2v1,ssig2w2,ssig2w1,ssig2w3,p11,p22,p12,p13,p21,p23,p31,p32,p33, F1,F2,F3,Q1,Q2,Q3,R1,R2,R3, beta_ll2,beta_ll1,beta_ll3,P_ll2,P_ll1,P_ll3,ppr1,pprl1_,pprl2_,pprl3_, lnL,itr,beta_tl11,beta_tl12,beta_tl13,beta_tl21,beta_tl22,beta_tl23,beta_tl31,beta_tl32,beta_tl33, P_tl11,P_tl12,P_tl13,P_tl21,P_tl22,P_tl23,P_tl31,P_tl32,P_tl33, pprl11,pprl12,pprl13,pprl21,pprl22,pprl23,pprl31,pprl32,pprl33, H_t2,H_t1,H_t3,eta_tl11,eta_tl12,eta_tl13,eta_tl21,eta_tl22,eta_tl23,eta_tl31,eta_tl32,eta_tl33, f_tl11,f_tl12,f_tl13,f_tl21,f_tl22,f_tl23,f_tl31,f_tl32,f_tl33, pdf11,pdf12,pdf13,pdf21,pdf22,pdf23,pdf31,pdf32,pdf33,pdf,lnf, Kt11,Kt12,Kt13,Kt21,Kt22,Kt23,Kt31,Kt32,Kt33, beta_tt11,beta_tt12,beta_tt13,beta_tt21,beta_tt22,beta_tt23,beta_tt31,beta_tt32,beta_tt33, P_tt11,P_tt12,P_tt13,P_tt21,P_tt22,P_tt23,P_tt31,P_tt32,P_tt33, beta_tt1, beta_tt2, beta_tt3,P_tt1,P_tt2,P_tt3,bb2,bb1,bb3,pprl, PR_TR,PPR_tl,J_pdf,gam,cor1,cor2,cor3,cov1,cov2,cov3,mu1,mu2,mu3; prmtr=trans(prmtr1); aa1=prmtr[1]; aa2=prmtr[2]; aa3=prmtr[3]; bb1=prmtr[4]; bb2=prmtr[5]; bb3=prmtr[6]; ssig2v1=prmtr[7]^2; ssig2v2=prmtr[8]^2; ssig2v3=prmtr[9]^2; ssig2w1=prmtr[10]^2; ssig2w2=prmtr[11]^2; ssig2w3=prmtr[12]^2; p11=prmtr[13];@Prb[ S(t)=1 | S(t-1)=1]@ p22=prmtr[14];@Prb[ S(t)=2 | S(t-1)=2]@ cor1=prmtr[15]; cor2=prmtr[16]; cor3=prmtr[17]; mu1=prmtr[18]; mu2=prmtr[19]; mu3=prmtr[20]; @ cor1=1; cor2=1; cor3=1; @ cov1=cor1*sqrt(ssig2v1*ssig2w1); cov2=cor2*sqrt(ssig2v2*ssig2w2); cov3=cor3*sqrt(ssig2v3*ssig2w3); @ F_j @ F1=(1~ 0 ~ 0 )| (0~aa1~bb1)| (0~ 1 ~ 0 ); F2=(1~ 0 ~ 0 )| (0~aa2~bb2)| (0~ 1 ~ 0 ); F3=(1~ 0 ~ 0 )| (0~aa3~bb3)| (0~ 1 ~ 0 ); @ Q_j @ Q1=( ssig2v1 ~ cov1 ~ 0 )| ( cov1 ~ ssig2w1 ~ 0 )| ( 0 ~ 0 ~ 0 ); Q2=( ssig2v2 ~ cov2 ~ 0 )| ( cov2 ~ ssig2w2 ~ 0 )| ( 0 ~ 0 ~ 0 ); Q3=( ssig2v3 ~ cov3 ~ 0 )| ( cov3 ~ ssig2w3 ~ 0 )| ( 0 ~ 0 ~ 0 ); R1=zeros(n,n); R2=zeros(n,n); R3=zeros(n,n); @ n by n @ beta_ll1=ymat[1]|0|0; beta_ll2=ymat[1]|0|0; beta_ll3=ymat[1]|0|0; p_ll1=500*eye(k); p_ll2=500*eye(k); p_ll3=500*eye(k); p_ll1=5*eye(k); p_ll2=5*eye(k); p_ll3=5*eye(k); p12=1-p11; p13=0; p21=0; p23=1-p22; p31=0; p32=0; p33=1; PR_TR= p11 ~ p12 ~ p13| p21 ~ p22 ~ p23| p31 ~ p32 ~ p33; PR_TR=vec(PR_TR); @ 9 by 1 @ @ At t=0, to calculate conditional probability St=1 , St=2 and St=3 @ pprl1_=0.99999; @ Unconditional Probability for St=1 @ pprl2_=0.00001; @ Unconditional Probability for St=2 @ pprl3_=0.00001; @ Unconditional Probability for St=3 @ pprl=pprl1_|pprl2_|pprl3_; @ 3 by 1 @ pprl=vec(pprl~pprl~pprl); @ 9 by 1 @ lnL=0; itr=1; do until itr>t; @ == Step 1 : Kalman Filter == @ @ beta_tl(i,j), S(t-1)=i, S(t)=j , k by 1 @ beta_tl11=mu1+ F1*beta_ll1; @ S(t-1)=1, S(t)=1 @ beta_tl21=mu1+ F1*beta_ll2; @ S(t-1)=2, S(t)=1 @ beta_tl31=mu1+ F1*beta_ll3; @ S(t-1)=3, S(t)=1 @ beta_tl12=mu2+ F2*beta_ll1; @ S(t-1)=1, S(t)=2 @ beta_tl22= mu2+F2*beta_ll2; @ S(t-1)=2, S(t)=2 @ beta_tl32=mu2+ F2*beta_ll3; @ S(t-1)=3, S(t)=2 @ beta_tl13= mu3+F3*beta_ll1; @ S(t-1)=1, S(t)=3 @ beta_tl23=mu3+ F3*beta_ll2; @ S(t-1)=2, S(t)=3 @ beta_tl33=mu3+ F3*beta_ll3; @ S(t-1)=3, S(t)=3 @ @ P_tl(i,j), S(t-1)=i, S(t)=j , k by k @ P_tl11=F1*P_ll1*F1' +Q1; P_tl21=F1*P_ll2*F1' +Q1; P_tl31=F1*P_ll3*F1' +Q1; P_tl12=F2*P_ll1*F2' +Q2; P_tl22=F2*P_ll2*F2' +Q2; P_tl32=F2*P_ll3*F2' +Q2; P_tl13=F3*P_ll1*F3' +Q3; P_tl23=F3*P_ll2*F3' +Q3; P_tl33=F3*P_ll3*F3' +Q3; @ H_t( j ), S(t)=j @ H_t1=1~1~0; @ 1 by k @ H_t2=1~1~0; H_t3=1~1~0; @ eta_tl(i,j), S(t-1)=i, S(t)=j , n by n @ eta_tl11=ymat[itr]-H_t1*beta_tl11 ; eta_tl21=ymat[itr]-H_t1*beta_tl21 ; eta_tl31=ymat[itr]-H_t1*beta_tl31 ; eta_tl12=ymat[itr]-H_t2*beta_tl12 ; eta_tl22=ymat[itr]-H_t2*beta_tl22 ; eta_tl32=ymat[itr]-H_t2*beta_tl32 ; eta_tl13=ymat[itr]-H_t3*beta_tl13 ; eta_tl23=ymat[itr]-H_t3*beta_tl23 ; eta_tl33=ymat[itr]-H_t3*beta_tl33 ; @ f_tl(i,j), S(t-1)=i, S(t)=j , n by n @ f_tl11=H_t1*P_tl11*H_t1'- R1; f_tl21=H_t1*P_tl21*H_t1'- R1; f_tl31=H_t1*P_tl31*H_t1'- R1; f_tl12=H_t2*P_tl12*H_t2'- R2; f_tl22=H_t2*P_tl22*H_t2'- R2; f_tl32=H_t2*P_tl32*H_t2'- R2; f_tl13=H_t3*P_tl13*H_t3'- R3; f_tl23=H_t3*P_tl23*H_t3'- R3; f_tl33=H_t3*P_tl33*H_t3'- R3; @ Kt(i,j), S(t-1)=i, S(t)=j , k by n @ Kt11=P_tl11*H_t1'*inv(f_tl11); Kt21=P_tl21*H_t1'*inv(f_tl21); Kt31=P_tl31*H_t1'*inv(f_tl31); Kt12=P_tl12*H_t2'*inv(f_tl12); Kt22=P_tl22*H_t2'*inv(f_tl22); Kt32=P_tl32*H_t2'*inv(f_tl32); Kt13=P_tl13*H_t3'*inv(f_tl13); Kt23=P_tl23*H_t3'*inv(f_tl23); Kt33=P_tl33*H_t3'*inv(f_tl33); @ == Updating == @ @ beta_tl(i,j), S(t-1)=i, S(t)=j , k by 1 @ beta_tt11= beta_tl11 + Kt11*eta_tl11; beta_tt21= beta_tl21 + Kt21*eta_tl21; beta_tt31= beta_tl31 + Kt31*eta_tl31; beta_tt12= beta_tl12 + Kt12*eta_tl12; beta_tt22= beta_tl22 + Kt22*eta_tl22; beta_tt32= beta_tl32 + Kt32*eta_tl32; beta_tt13= beta_tl13 + Kt13*eta_tl13; beta_tt23= beta_tl23 + Kt23*eta_tl23; beta_tt33= beta_tl33 + Kt33*eta_tl33; @ P_tl(i,j), S(t-1)=i, S(t)=j , k by k @ P_tt11=(eye(k)-P_tl11*H_t1'*inv(f_tl11)*H_t1)*P_tl11; P_tt21=(eye(k)-P_tl21*H_t1'*inv(f_tl21)*H_t1)*P_tl21; P_tt31=(eye(k)-P_tl31*H_t1'*inv(f_tl31)*H_t1)*P_tl31; P_tt12=(eye(k)-P_tl12*H_t2'*inv(f_tl12)*H_t2)*P_tl12; P_tt22=(eye(k)-P_tl22*H_t2'*inv(f_tl22)*H_t2)*P_tl22; P_tt32=(eye(k)-P_tl32*H_t2'*inv(f_tl32)*H_t2)*P_tl32; P_tt13=(eye(k)-P_tl13*H_t3'*inv(f_tl13)*H_t3)*P_tl13; P_tt23=(eye(k)-P_tl23*H_t3'*inv(f_tl23)*H_t3)*P_tl23; P_tt33=(eye(k)-P_tl33*H_t3'*inv(f_tl33)*H_t3)*P_tl33; @ == Step 2 : Hamilton Filter == @ PPR_tl=pprl.* PR_TR; @ 9 by 1 @ pdf11=(1/sqrt(2*pi*det(f_tl11)))*exp(-0.5*eta_tl11'inv(f_tl11)*eta_tl11); pdf21=(1/sqrt(2*pi*det(f_tl21)))*exp(-0.5*eta_tl21'inv(f_tl21)*eta_tl21); pdf31=(1/sqrt(2*pi*det(f_tl31)))*exp(-0.5*eta_tl31'inv(f_tl31)*eta_tl31); pdf12=(1/sqrt(2*pi*det(f_tl12)))*exp(-0.5*eta_tl12'inv(f_tl12)*eta_tl12); pdf22=(1/sqrt(2*pi*det(f_tl22)))*exp(-0.5*eta_tl22'inv(f_tl22)*eta_tl22); pdf32=(1/sqrt(2*pi*det(f_tl32)))*exp(-0.5*eta_tl32'inv(f_tl32)*eta_tl32); pdf13=(1/sqrt(2*pi*det(f_tl13)))*exp(-0.5*eta_tl13'inv(f_tl13)*eta_tl13); pdf23=(1/sqrt(2*pi*det(f_tl23)))*exp(-0.5*eta_tl23'inv(f_tl23)*eta_tl23); pdf33=(1/sqrt(2*pi*det(f_tl33)))*exp(-0.5*eta_tl33'inv(f_tl33)*eta_tl33); pdf=pdf11|pdf21|pdf31|pdf12|pdf22|pdf32|pdf13|pdf23|pdf33; pdf=pdf.*PPR_tl; J_pdf=sumc(pdf); if itr>start; lnf=ln(J_pdf); else; lnf=0; endif; lnL=lnL+lnf; pprl1_=sumc(pdf[1:3])/J_pdf; pprl2_=sumc(pdf[4:6])/J_pdf; pprl3_=sumc(pdf[7:9])/J_pdf; @ == Step 3 : Approximations == @ @ beta_tt(j), S(t) = j, k by 1 @ beta_tt1=(beta_tt11*(pdf[1]/j_pdf)+beta_tt21*(pdf[2]/j_pdf)+beta_tt31*(pdf[3]/j_pdf))/pprl1_; @ beta_tt(1) @ beta_tt2=(beta_tt12*(pdf[4]/j_pdf)+beta_tt22*(pdf[5]/j_pdf)+beta_tt32*(pdf[6]/j_pdf))/pprl2_; @ beta_tt(2) @ beta_tt3=(beta_tt13*(pdf[7]/j_pdf)+beta_tt23*(pdf[8]/j_pdf)+beta_tt33*(pdf[9]/j_pdf))/pprl3_; @ beta_tt(3) @ @ P_tt(j), S(t) = j, k by k @ P_tt1=( (pdf[1]/J_pdf)*(P_tt11+(beta_tt1-beta_tt11)*(beta_tt1-beta_tt11)') +(pdf[2]/J_pdf)*(P_tt21+(beta_tt1-beta_tt21)*(beta_tt1-beta_tt21)') +(pdf[3]/J_pdf)*(P_tt31+(beta_tt1-beta_tt31)*(beta_tt1-beta_tt31)') )/pprl1_; @ P_tt(1) @ P_tt2=( (pdf[4]/J_pdf)*(P_tt12+(beta_tt2-beta_tt12)*(beta_tt2-beta_tt12)') +(pdf[5]/J_pdf)*(P_tt22+(beta_tt2-beta_tt22)*(beta_tt2-beta_tt22)') +(pdf[6]/J_pdf)*(P_tt32+(beta_tt2-beta_tt32)*(beta_tt2-beta_tt32)') )/pprl2_; @ P_tt(2) @ P_tt3=( (pdf[7]/J_pdf)*(P_tt13+(beta_tt3-beta_tt13)*(beta_tt3-beta_tt13)') +(pdf[8]/J_pdf)*(P_tt23+(beta_tt3-beta_tt23)*(beta_tt3-beta_tt23)') +(pdf[9]/J_pdf)*(P_tt33+(beta_tt3-beta_tt33)*(beta_tt3-beta_tt33)') )/pprl3_; @ P_tt(3) @ pprl=pprl1_|pprl2_|pprl3_; @ 3 by 1 @ pprl=vec(pprl~pprl~pprl); @ 9 by 1 @ beta_ll1=beta_tt1; beta_ll2= beta_tt2; beta_ll3= beta_tt3; P_ll1= P_tt1;P_ll2= P_tt2;P_ll3= P_tt3; itr=itr+1; endo; retp(-lnL); endp; /*================= Lik_f0===================*/ proc lik_f0(prmtr1); local prmtr,mmu0,mmu1,aa3,aa1,aa2,ssig2v2,ssig2v3,ssig2v1,ssig2w2,ssig2w1,ssig2w3,p11,p22,p12,p13,p21,p23,p31,p32,p33, F1,F2,F3,Q1,Q2,Q3,R1,R2,R3, beta_ll2,beta_ll1,beta_ll3,P_ll2,P_ll1,P_ll3,ppr1,pprl1_,pprl2_,pprl3_, lnL,itr,beta_tl11,beta_tl12,beta_tl13,beta_tl21,beta_tl22,beta_tl23,beta_tl31,beta_tl32,beta_tl33, P_tl11,P_tl12,P_tl13,P_tl21,P_tl22,P_tl23,P_tl31,P_tl32,P_tl33, pprl11,pprl12,pprl13,pprl21,pprl22,pprl23,pprl31,pprl32,pprl33, H_t2,H_t1,H_t3,eta_tl11,eta_tl12,eta_tl13,eta_tl21,eta_tl22,eta_tl23,eta_tl31,eta_tl32,eta_tl33, f_tl11,f_tl12,f_tl13,f_tl21,f_tl22,f_tl23,f_tl31,f_tl32,f_tl33, pdf11,pdf12,pdf13,pdf21,pdf22,pdf23,pdf31,pdf32,pdf33,pdf,lnf, Kt11,Kt12,Kt13,Kt21,Kt22,Kt23,Kt31,Kt32,Kt33, beta_tt11,beta_tt12,beta_tt13,beta_tt21,beta_tt22,beta_tt23,beta_tt31,beta_tt32,beta_tt33, P_tt11,P_tt12,P_tt13,P_tt21,P_tt22,P_tt23,P_tt31,P_tt32,P_tt33, beta_tt1, beta_tt2, beta_tt3,P_tt1,P_tt2,P_tt3,bb2,bb1,bb3,pprl, PR_TR,PPR_tl,J_pdf,gam,cor1,cor2,cor3,cov1,cov2,cov3,mu1,mu2,mu3; prmtr=trans0(prmtr1); aa1=prmtr[1]; aa2=prmtr[2]; aa3=prmtr[3]; bb1=prmtr[4]; bb2=prmtr[5]; bb3=prmtr[6]; ssig2v1=prmtr[7]^2; ssig2v2=prmtr[8]^2; ssig2v3=prmtr[9]^2; ssig2w1=prmtr[10]^2; ssig2w2=ssig2w1; ssig2w3=ssig2w1; p11=prmtr[11];@Prb[ S(t)=1 | S(t-1)=1]@ p22=prmtr[12];@Prb[ S(t)=0 | S(t-1)=0]@ cor1=prmtr[13]; cor2=prmtr[14]; cor3=prmtr[15]; mu1=prmtr[16]; mu2=prmtr[17]; mu3=prmtr[18]; cov1=cor1*sqrt(ssig2v1*ssig2w1); cov2=cor2*sqrt(ssig2v2*ssig2w2); cov3=cor3*sqrt(ssig2v3*ssig2w3); @ F_j @ F1=(1~ 0 ~ 0 )| (0~aa1~bb1)| (0~ 1 ~ 0 ); F2=(1~ 0 ~ 0 )| (0~aa2~bb2)| (0~ 1 ~ 0 ); F3=(1~ 0 ~ 0 )| (0~aa3~bb3)| (0~ 1 ~ 0 ); @ Q_j @ Q1=( ssig2v1 ~ cov1 ~ 0 )| ( cov1 ~ ssig2w1 ~ 0 )| ( 0 ~ 0 ~ 0 ); Q2=( ssig2v2 ~ cov2 ~ 0 )| ( cov2 ~ ssig2w2 ~ 0 )| ( 0 ~ 0 ~ 0 ); Q3=( ssig2v3 ~ cov3 ~ 0 )| ( cov3 ~ ssig2w3 ~ 0 )| ( 0 ~ 0 ~ 0 ); R1=zeros(n,n); R2=zeros(n,n); R3=zeros(n,n); @ n by n @ beta_ll1=ymat[1]|0|0; beta_ll2=ymat[1]|0|0; beta_ll3=ymat[1]|0|0; p_ll1=500*eye(k); p_ll2=500*eye(k); p_ll3=500*eye(k); p_ll1=5*eye(k); p_ll2=5*eye(k); p_ll3=5*eye(k); p12=1-p11; p13=0; p21=0; p23=1-p22; p31=0; p32=0; p33=1; PR_TR= p11 ~ p12 ~ p13| p21 ~ p22 ~ p23| p31 ~ p32 ~ p33; PR_TR=vec(PR_TR); @ 9 by 1 @ @ At t=0, to calculate conditional probability St=1 , St=2 and St=3 @ pprl1_=0.99999; @ Unconditional Probability for St=1 @ pprl2_=0.00001; @ Unconditional Probability for St=2 @ pprl3_=0.00001; @ Unconditional Probability for St=3 @ pprl=pprl1_|pprl2_|pprl3_; @ 3 by 1 @ pprl=vec(pprl~pprl~pprl); @ 9 by 1 @ lnL=0; itr=1; do until itr>t; @ == Step 1 : Kalman Filter == @ @ beta_tl(i,j), S(t-1)=i, S(t)=j , k by 1 @ beta_tl11=mu1+ F1*beta_ll1; @ S(t-1)=1, S(t)=1 @ beta_tl21=mu1+ F1*beta_ll2; @ S(t-1)=2, S(t)=1 @ beta_tl31=mu1+ F1*beta_ll3; @ S(t-1)=3, S(t)=1 @ beta_tl12=mu2+ F2*beta_ll1; @ S(t-1)=1, S(t)=2 @ beta_tl22= mu2+F2*beta_ll2; @ S(t-1)=2, S(t)=2 @ beta_tl32=mu2+ F2*beta_ll3; @ S(t-1)=3, S(t)=2 @ beta_tl13= mu3+F3*beta_ll1; @ S(t-1)=1, S(t)=3 @ beta_tl23=mu3+ F3*beta_ll2; @ S(t-1)=2, S(t)=3 @ beta_tl33=mu3+ F3*beta_ll3; @ S(t-1)=3, S(t)=3 @ @ P_tl(i,j), S(t-1)=i, S(t)=j , k by k @ P_tl11=F1*P_ll1*F1' +Q1; P_tl21=F1*P_ll2*F1' +Q1; P_tl31=F1*P_ll3*F1' +Q1; P_tl12=F2*P_ll1*F2' +Q2; P_tl22=F2*P_ll2*F2' +Q2; P_tl32=F2*P_ll3*F2' +Q2; P_tl13=F3*P_ll1*F3' +Q3; P_tl23=F3*P_ll2*F3' +Q3; P_tl33=F3*P_ll3*F3' +Q3; @ H_t( j ), S(t)=j @ H_t1=1~1~0; @ 1 by k @ H_t2=1~1~0; H_t3=1~1~0; @ eta_tl(i,j), S(t-1)=i, S(t)=j , n by n @ eta_tl11=ymat[itr]-H_t1*beta_tl11 ; eta_tl21=ymat[itr]-H_t1*beta_tl21 ; eta_tl31=ymat[itr]-H_t1*beta_tl31 ; eta_tl12=ymat[itr]-H_t2*beta_tl12 ; eta_tl22=ymat[itr]-H_t2*beta_tl22 ; eta_tl32=ymat[itr]-H_t2*beta_tl32 ; eta_tl13=ymat[itr]-H_t3*beta_tl13 ; eta_tl23=ymat[itr]-H_t3*beta_tl23 ; eta_tl33=ymat[itr]-H_t3*beta_tl33 ; @ f_tl(i,j), S(t-1)=i, S(t)=j , n by n @ f_tl11=H_t1*P_tl11*H_t1'- R1; f_tl21=H_t1*P_tl21*H_t1'- R1; f_tl31=H_t1*P_tl31*H_t1'- R1; f_tl12=H_t2*P_tl12*H_t2'- R2; f_tl22=H_t2*P_tl22*H_t2'- R2; f_tl32=H_t2*P_tl32*H_t2'- R2; f_tl13=H_t3*P_tl13*H_t3'- R3; f_tl23=H_t3*P_tl23*H_t3'- R3; f_tl33=H_t3*P_tl33*H_t3'- R3; @ Kt(i,j), S(t-1)=i, S(t)=j , k by n @ Kt11=P_tl11*H_t1'*inv(f_tl11); Kt21=P_tl21*H_t1'*inv(f_tl21); Kt31=P_tl31*H_t1'*inv(f_tl31); Kt12=P_tl12*H_t2'*inv(f_tl12); Kt22=P_tl22*H_t2'*inv(f_tl22); Kt32=P_tl32*H_t2'*inv(f_tl32); Kt13=P_tl13*H_t3'*inv(f_tl13); Kt23=P_tl23*H_t3'*inv(f_tl23); Kt33=P_tl33*H_t3'*inv(f_tl33); @ == Updating == @ @ beta_tl(i,j), S(t-1)=i, S(t)=j , k by 1 @ beta_tt11= beta_tl11 + Kt11*eta_tl11; beta_tt21= beta_tl21 + Kt21*eta_tl21; beta_tt31= beta_tl31 + Kt31*eta_tl31; beta_tt12= beta_tl12 + Kt12*eta_tl12; beta_tt22= beta_tl22 + Kt22*eta_tl22; beta_tt32= beta_tl32 + Kt32*eta_tl32; beta_tt13= beta_tl13 + Kt13*eta_tl13; beta_tt23= beta_tl23 + Kt23*eta_tl23; beta_tt33= beta_tl33 + Kt33*eta_tl33; @ P_tl(i,j), S(t-1)=i, S(t)=j , k by k @ P_tt11=(eye(k)-P_tl11*H_t1'*inv(f_tl11)*H_t1)*P_tl11; P_tt21=(eye(k)-P_tl21*H_t1'*inv(f_tl21)*H_t1)*P_tl21; P_tt31=(eye(k)-P_tl31*H_t1'*inv(f_tl31)*H_t1)*P_tl31; P_tt12=(eye(k)-P_tl12*H_t2'*inv(f_tl12)*H_t2)*P_tl12; P_tt22=(eye(k)-P_tl22*H_t2'*inv(f_tl22)*H_t2)*P_tl22; P_tt32=(eye(k)-P_tl32*H_t2'*inv(f_tl32)*H_t2)*P_tl32; P_tt13=(eye(k)-P_tl13*H_t3'*inv(f_tl13)*H_t3)*P_tl13; P_tt23=(eye(k)-P_tl23*H_t3'*inv(f_tl23)*H_t3)*P_tl23; P_tt33=(eye(k)-P_tl33*H_t3'*inv(f_tl33)*H_t3)*P_tl33; @ == Step 2 : Hamilton Filter == @ PPR_tl=pprl.* PR_TR; @ 9 by 1 @ pdf11=(1/sqrt(2*pi*det(f_tl11)))*exp(-0.5*eta_tl11'inv(f_tl11)*eta_tl11); pdf21=(1/sqrt(2*pi*det(f_tl21)))*exp(-0.5*eta_tl21'inv(f_tl21)*eta_tl21); pdf31=(1/sqrt(2*pi*det(f_tl31)))*exp(-0.5*eta_tl31'inv(f_tl31)*eta_tl31); pdf12=(1/sqrt(2*pi*det(f_tl12)))*exp(-0.5*eta_tl12'inv(f_tl12)*eta_tl12); pdf22=(1/sqrt(2*pi*det(f_tl22)))*exp(-0.5*eta_tl22'inv(f_tl22)*eta_tl22); pdf32=(1/sqrt(2*pi*det(f_tl32)))*exp(-0.5*eta_tl32'inv(f_tl32)*eta_tl32); pdf13=(1/sqrt(2*pi*det(f_tl13)))*exp(-0.5*eta_tl13'inv(f_tl13)*eta_tl13); pdf23=(1/sqrt(2*pi*det(f_tl23)))*exp(-0.5*eta_tl23'inv(f_tl23)*eta_tl23); pdf33=(1/sqrt(2*pi*det(f_tl33)))*exp(-0.5*eta_tl33'inv(f_tl33)*eta_tl33); pdf=pdf11|pdf21|pdf31|pdf12|pdf22|pdf32|pdf13|pdf23|pdf33; pdf=pdf.*PPR_tl; J_pdf=sumc(pdf); if itr>start; lnf=ln(J_pdf); else; lnf=0; endif; lnL=lnL+lnf; pprl1_=sumc(pdf[1:3])/J_pdf; pprl2_=sumc(pdf[4:6])/J_pdf; pprl3_=sumc(pdf[7:9])/J_pdf; @ == Step 3 : Approximations == @ @ beta_tt(j), S(t) = j, k by 1 @ beta_tt1=(beta_tt11*(pdf[1]/j_pdf)+beta_tt21*(pdf[2]/j_pdf)+beta_tt31*(pdf[3]/j_pdf))/pprl1_; @ beta_tt(1) @ beta_tt2=(beta_tt12*(pdf[4]/j_pdf)+beta_tt22*(pdf[5]/j_pdf)+beta_tt32*(pdf[6]/j_pdf))/pprl2_; @ beta_tt(2) @ beta_tt3=(beta_tt13*(pdf[7]/j_pdf)+beta_tt23*(pdf[8]/j_pdf)+beta_tt33*(pdf[9]/j_pdf))/pprl3_; @ beta_tt(3) @ @ P_tt(j), S(t) = j, k by k @ P_tt1=( (pdf[1]/J_pdf)*(P_tt11+(beta_tt1-beta_tt11)*(beta_tt1-beta_tt11)') +(pdf[2]/J_pdf)*(P_tt21+(beta_tt1-beta_tt21)*(beta_tt1-beta_tt21)') +(pdf[3]/J_pdf)*(P_tt31+(beta_tt1-beta_tt31)*(beta_tt1-beta_tt31)') )/pprl1_; @ P_tt(1) @ P_tt2=( (pdf[4]/J_pdf)*(P_tt12+(beta_tt2-beta_tt12)*(beta_tt2-beta_tt12)') +(pdf[5]/J_pdf)*(P_tt22+(beta_tt2-beta_tt22)*(beta_tt2-beta_tt22)') +(pdf[6]/J_pdf)*(P_tt32+(beta_tt2-beta_tt32)*(beta_tt2-beta_tt32)') )/pprl2_; @ P_tt(2) @ P_tt3=( (pdf[7]/J_pdf)*(P_tt13+(beta_tt3-beta_tt13)*(beta_tt3-beta_tt13)') +(pdf[8]/J_pdf)*(P_tt23+(beta_tt3-beta_tt23)*(beta_tt3-beta_tt23)') +(pdf[9]/J_pdf)*(P_tt33+(beta_tt3-beta_tt33)*(beta_tt3-beta_tt33)') )/pprl3_; @ P_tt(3) @ pprl=pprl1_|pprl2_|pprl3_; @ 3 by 1 @ pprl=vec(pprl~pprl~pprl); @ 9 by 1 @ beta_ll1=beta_tt1; beta_ll2= beta_tt2; beta_ll3= beta_tt3; P_ll1= P_tt1;P_ll2= P_tt2;P_ll3= P_tt3; itr=itr+1; endo; retp(-lnL); endp; /*================= Lik_f0===================*/ proc lik_f1(prmtr1); local prmtr,mmu0,mmu1,aa3,aa1,aa2,ssig2v2,ssig2v3,ssig2v1,ssig2w2,ssig2w1,ssig2w3,p11,p22,p12,p13,p21,p23,p31,p32,p33, F1,F2,F3,Q1,Q2,Q3,R1,R2,R3, beta_ll2,beta_ll1,beta_ll3,P_ll2,P_ll1,P_ll3,ppr1,pprl1_,pprl2_,pprl3_, lnL,itr,beta_tl11,beta_tl12,beta_tl13,beta_tl21,beta_tl22,beta_tl23,beta_tl31,beta_tl32,beta_tl33, P_tl11,P_tl12,P_tl13,P_tl21,P_tl22,P_tl23,P_tl31,P_tl32,P_tl33, pprl11,pprl12,pprl13,pprl21,pprl22,pprl23,pprl31,pprl32,pprl33, H_t2,H_t1,H_t3,eta_tl11,eta_tl12,eta_tl13,eta_tl21,eta_tl22,eta_tl23,eta_tl31,eta_tl32,eta_tl33, f_tl11,f_tl12,f_tl13,f_tl21,f_tl22,f_tl23,f_tl31,f_tl32,f_tl33, pdf11,pdf12,pdf13,pdf21,pdf22,pdf23,pdf31,pdf32,pdf33,pdf,lnf, Kt11,Kt12,Kt13,Kt21,Kt22,Kt23,Kt31,Kt32,Kt33, beta_tt11,beta_tt12,beta_tt13,beta_tt21,beta_tt22,beta_tt23,beta_tt31,beta_tt32,beta_tt33, P_tt11,P_tt12,P_tt13,P_tt21,P_tt22,P_tt23,P_tt31,P_tt32,P_tt33, beta_tt1, beta_tt2, beta_tt3,P_tt1,P_tt2,P_tt3,bb2,bb1,bb3,pprl, PR_TR,PPR_tl,J_pdf,gam,cor1,cor2,cor3,cov1,cov2,cov3,mu1,mu2,mu3; prmtr=trans0(prmtr1); aa1=prmtr[1]; aa2=prmtr[2]; aa3=prmtr[3]; bb1=prmtr[4]; bb2=prmtr[5]; bb3=prmtr[6]; ssig2v1=prmtr[7]^2; ssig2v2=ssig2v1; ssig2v3=ssig2v1; ssig2w1=prmtr[8]^2; ssig2w2=prmtr[9]^2; ssig2w3=prmtr[10]^2; p11=prmtr[11];@Prb[ S(t)=1 | S(t-1)=1]@ p22=prmtr[12];@Prb[ S(t)=0 | S(t-1)=0]@ cor1=prmtr[13]; cor2=prmtr[14]; cor3=prmtr[15]; mu1=prmtr[16]; mu2=prmtr[17]; mu3=prmtr[18]; cov1=cor1*sqrt(ssig2v1*ssig2w1); cov2=cor2*sqrt(ssig2v2*ssig2w2); cov3=cor3*sqrt(ssig2v3*ssig2w3); @ F_j @ F1=(1~ 0 ~ 0 )| (0~aa1~bb1)| (0~ 1 ~ 0 ); F2=(1~ 0 ~ 0 )| (0~aa2~bb2)| (0~ 1 ~ 0 ); F3=(1~ 0 ~ 0 )| (0~aa3~bb3)| (0~ 1 ~ 0 ); @ Q_j @ Q1=( ssig2v1 ~ cov1 ~ 0 )| ( cov1 ~ ssig2w1 ~ 0 )| ( 0 ~ 0 ~ 0 ); Q2=( ssig2v2 ~ cov2 ~ 0 )| ( cov2 ~ ssig2w2 ~ 0 )| ( 0 ~ 0 ~ 0 ); Q3=( ssig2v3 ~ cov3 ~ 0 )| ( cov3 ~ ssig2w3 ~ 0 )| ( 0 ~ 0 ~ 0 ); R1=zeros(n,n); R2=zeros(n,n); R3=zeros(n,n); @ n by n @ beta_ll1=ymat[1]|0|0; beta_ll2=ymat[1]|0|0; beta_ll3=ymat[1]|0|0; p_ll1=500*eye(k); p_ll2=500*eye(k); p_ll3=500*eye(k); p_ll1=5*eye(k); p_ll2=5*eye(k); p_ll3=5*eye(k); p12=1-p11; p13=0; p21=0; p23=1-p22; p31=0; p32=0; p33=1; PR_TR= p11 ~ p12 ~ p13| p21 ~ p22 ~ p23| p31 ~ p32 ~ p33; PR_TR=vec(PR_TR); @ 9 by 1 @ @ At t=0, to calculate conditional probability St=1 , St=2 and St=3 @ pprl1_=0.99999; @ Unconditional Probability for St=1 @ pprl2_=0.00001; @ Unconditional Probability for St=2 @ pprl3_=0.00001; @ Unconditional Probability for St=3 @ pprl=pprl1_|pprl2_|pprl3_; @ 3 by 1 @ pprl=vec(pprl~pprl~pprl); @ 9 by 1 @ lnL=0; itr=1; do until itr>t; @ == Step 1 : Kalman Filter == @ @ beta_tl(i,j), S(t-1)=i, S(t)=j , k by 1 @ beta_tl11=mu1+ F1*beta_ll1; @ S(t-1)=1, S(t)=1 @ beta_tl21=mu1+ F1*beta_ll2; @ S(t-1)=2, S(t)=1 @ beta_tl31=mu1+ F1*beta_ll3; @ S(t-1)=3, S(t)=1 @ beta_tl12=mu2+ F2*beta_ll1; @ S(t-1)=1, S(t)=2 @ beta_tl22= mu2+F2*beta_ll2; @ S(t-1)=2, S(t)=2 @ beta_tl32=mu2+ F2*beta_ll3; @ S(t-1)=3, S(t)=2 @ beta_tl13= mu3+F3*beta_ll1; @ S(t-1)=1, S(t)=3 @ beta_tl23=mu3+ F3*beta_ll2; @ S(t-1)=2, S(t)=3 @ beta_tl33=mu3+ F3*beta_ll3; @ S(t-1)=3, S(t)=3 @ @ P_tl(i,j), S(t-1)=i, S(t)=j , k by k @ P_tl11=F1*P_ll1*F1' +Q1; P_tl21=F1*P_ll2*F1' +Q1; P_tl31=F1*P_ll3*F1' +Q1; P_tl12=F2*P_ll1*F2' +Q2; P_tl22=F2*P_ll2*F2' +Q2; P_tl32=F2*P_ll3*F2' +Q2; P_tl13=F3*P_ll1*F3' +Q3; P_tl23=F3*P_ll2*F3' +Q3; P_tl33=F3*P_ll3*F3' +Q3; @ H_t( j ), S(t)=j @ H_t1=1~1~0; @ 1 by k @ H_t2=1~1~0; H_t3=1~1~0; @ eta_tl(i,j), S(t-1)=i, S(t)=j , n by n @ eta_tl11=ymat[itr]-H_t1*beta_tl11 ; eta_tl21=ymat[itr]-H_t1*beta_tl21 ; eta_tl31=ymat[itr]-H_t1*beta_tl31 ; eta_tl12=ymat[itr]-H_t2*beta_tl12 ; eta_tl22=ymat[itr]-H_t2*beta_tl22 ; eta_tl32=ymat[itr]-H_t2*beta_tl32 ; eta_tl13=ymat[itr]-H_t3*beta_tl13 ; eta_tl23=ymat[itr]-H_t3*beta_tl23 ; eta_tl33=ymat[itr]-H_t3*beta_tl33 ; @ f_tl(i,j), S(t-1)=i, S(t)=j , n by n @ f_tl11=H_t1*P_tl11*H_t1'- R1; f_tl21=H_t1*P_tl21*H_t1'- R1; f_tl31=H_t1*P_tl31*H_t1'- R1; f_tl12=H_t2*P_tl12*H_t2'- R2; f_tl22=H_t2*P_tl22*H_t2'- R2; f_tl32=H_t2*P_tl32*H_t2'- R2; f_tl13=H_t3*P_tl13*H_t3'- R3; f_tl23=H_t3*P_tl23*H_t3'- R3; f_tl33=H_t3*P_tl33*H_t3'- R3; @ Kt(i,j), S(t-1)=i, S(t)=j , k by n @ Kt11=P_tl11*H_t1'*inv(f_tl11); Kt21=P_tl21*H_t1'*inv(f_tl21); Kt31=P_tl31*H_t1'*inv(f_tl31); Kt12=P_tl12*H_t2'*inv(f_tl12); Kt22=P_tl22*H_t2'*inv(f_tl22); Kt32=P_tl32*H_t2'*inv(f_tl32); Kt13=P_tl13*H_t3'*inv(f_tl13); Kt23=P_tl23*H_t3'*inv(f_tl23); Kt33=P_tl33*H_t3'*inv(f_tl33); @ == Updating == @ @ beta_tl(i,j), S(t-1)=i, S(t)=j , k by 1 @ beta_tt11= beta_tl11 + Kt11*eta_tl11; beta_tt21= beta_tl21 + Kt21*eta_tl21; beta_tt31= beta_tl31 + Kt31*eta_tl31; beta_tt12= beta_tl12 + Kt12*eta_tl12; beta_tt22= beta_tl22 + Kt22*eta_tl22; beta_tt32= beta_tl32 + Kt32*eta_tl32; beta_tt13= beta_tl13 + Kt13*eta_tl13; beta_tt23= beta_tl23 + Kt23*eta_tl23; beta_tt33= beta_tl33 + Kt33*eta_tl33; @ P_tl(i,j), S(t-1)=i, S(t)=j , k by k @ P_tt11=(eye(k)-P_tl11*H_t1'*inv(f_tl11)*H_t1)*P_tl11; P_tt21=(eye(k)-P_tl21*H_t1'*inv(f_tl21)*H_t1)*P_tl21; P_tt31=(eye(k)-P_tl31*H_t1'*inv(f_tl31)*H_t1)*P_tl31; P_tt12=(eye(k)-P_tl12*H_t2'*inv(f_tl12)*H_t2)*P_tl12; P_tt22=(eye(k)-P_tl22*H_t2'*inv(f_tl22)*H_t2)*P_tl22; P_tt32=(eye(k)-P_tl32*H_t2'*inv(f_tl32)*H_t2)*P_tl32; P_tt13=(eye(k)-P_tl13*H_t3'*inv(f_tl13)*H_t3)*P_tl13; P_tt23=(eye(k)-P_tl23*H_t3'*inv(f_tl23)*H_t3)*P_tl23; P_tt33=(eye(k)-P_tl33*H_t3'*inv(f_tl33)*H_t3)*P_tl33; @ == Step 2 : Hamilton Filter == @ PPR_tl=pprl.* PR_TR; @ 9 by 1 @ pdf11=(1/sqrt(2*pi*det(f_tl11)))*exp(-0.5*eta_tl11'inv(f_tl11)*eta_tl11); pdf21=(1/sqrt(2*pi*det(f_tl21)))*exp(-0.5*eta_tl21'inv(f_tl21)*eta_tl21); pdf31=(1/sqrt(2*pi*det(f_tl31)))*exp(-0.5*eta_tl31'inv(f_tl31)*eta_tl31); pdf12=(1/sqrt(2*pi*det(f_tl12)))*exp(-0.5*eta_tl12'inv(f_tl12)*eta_tl12); pdf22=(1/sqrt(2*pi*det(f_tl22)))*exp(-0.5*eta_tl22'inv(f_tl22)*eta_tl22); pdf32=(1/sqrt(2*pi*det(f_tl32)))*exp(-0.5*eta_tl32'inv(f_tl32)*eta_tl32); pdf13=(1/sqrt(2*pi*det(f_tl13)))*exp(-0.5*eta_tl13'inv(f_tl13)*eta_tl13); pdf23=(1/sqrt(2*pi*det(f_tl23)))*exp(-0.5*eta_tl23'inv(f_tl23)*eta_tl23); pdf33=(1/sqrt(2*pi*det(f_tl33)))*exp(-0.5*eta_tl33'inv(f_tl33)*eta_tl33); pdf=pdf11|pdf21|pdf31|pdf12|pdf22|pdf32|pdf13|pdf23|pdf33; pdf=pdf.*PPR_tl; J_pdf=sumc(pdf); if itr>start; lnf=ln(J_pdf); else; lnf=0; endif; lnL=lnL+lnf; pprl1_=sumc(pdf[1:3])/J_pdf; pprl2_=sumc(pdf[4:6])/J_pdf; pprl3_=sumc(pdf[7:9])/J_pdf; @ == Step 3 : Approximations == @ @ beta_tt(j), S(t) = j, k by 1 @ beta_tt1=(beta_tt11*(pdf[1]/j_pdf)+beta_tt21*(pdf[2]/j_pdf)+beta_tt31*(pdf[3]/j_pdf))/pprl1_; @ beta_tt(1) @ beta_tt2=(beta_tt12*(pdf[4]/j_pdf)+beta_tt22*(pdf[5]/j_pdf)+beta_tt32*(pdf[6]/j_pdf))/pprl2_; @ beta_tt(1) @ beta_tt3=(beta_tt13*(pdf[7]/j_pdf)+beta_tt23*(pdf[8]/j_pdf)+beta_tt33*(pdf[9]/j_pdf))/pprl3_; @ beta_tt(1) @ @ P_tt(j), S(t) = j, k by k @ P_tt1=( (pdf[1]/J_pdf)*(P_tt11+(beta_tt1-beta_tt11)*(beta_tt1-beta_tt11)') +(pdf[2]/J_pdf)*(P_tt21+(beta_tt1-beta_tt21)*(beta_tt1-beta_tt21)') +(pdf[3]/J_pdf)*(P_tt31+(beta_tt1-beta_tt31)*(beta_tt1-beta_tt31)') )/pprl1_; @ P_tt(1) @ P_tt2=( (pdf[4]/J_pdf)*(P_tt12+(beta_tt2-beta_tt12)*(beta_tt2-beta_tt12)') +(pdf[5]/J_pdf)*(P_tt22+(beta_tt2-beta_tt22)*(beta_tt2-beta_tt22)') +(pdf[6]/J_pdf)*(P_tt32+(beta_tt2-beta_tt32)*(beta_tt2-beta_tt32)') )/pprl2_; @ P_tt(1) @ P_tt3=( (pdf[7]/J_pdf)*(P_tt13+(beta_tt3-beta_tt13)*(beta_tt3-beta_tt13)') +(pdf[8]/J_pdf)*(P_tt23+(beta_tt3-beta_tt23)*(beta_tt3-beta_tt23)') +(pdf[9]/J_pdf)*(P_tt33+(beta_tt3-beta_tt33)*(beta_tt3-beta_tt33)') )/pprl3_; @ P_tt(1) @ pprl=pprl1_|pprl2_|pprl3_; @ 3 by 1 @ pprl=vec(pprl~pprl~pprl); @ 9 by 1 @ beta_ll1=beta_tt1; beta_ll2= beta_tt2; beta_ll3= beta_tt3; P_ll1= P_tt1;P_ll2= P_tt2;P_ll3= P_tt3; itr=itr+1; endo; retp(-lnL); endp; /*================= Filtering ===================*/ proc(7)=filter(prmtr1); local prmtr,mmu0,mmu1,aa3,aa1,aa2,ssig2v2,ssig2v3,ssig2v1,ssig2w2,ssig2w1,ssig2w3,p11,p22,p12,p13,p21,p23,p31,p32,p33, F1,F2,F3,Q1,Q2,Q3,R1,R2,R3, beta_ll2,beta_ll1,beta_ll3,P_ll2,P_ll1,P_ll3,ppr1,pprl1_,pprl2_,pprl3_, lnL,itr,beta_tl11,beta_tl12,beta_tl13,beta_tl21,beta_tl22,beta_tl23,beta_tl31,beta_tl32,beta_tl33, P_tl11,P_tl12,P_tl13,P_tl21,P_tl22,P_tl23,P_tl31,P_tl32,P_tl33, pprl11,pprl12,pprl13,pprl21,pprl22,pprl23,pprl31,pprl32,pprl33, H_t2,H_t1,H_t3,eta_tl11,eta_tl12,eta_tl13,eta_tl21,eta_tl22,eta_tl23,eta_tl31,eta_tl32,eta_tl33, f_tl11,f_tl12,f_tl13,f_tl21,f_tl22,f_tl23,f_tl31,f_tl32,f_tl33, pdf11,pdf12,pdf13,pdf21,pdf22,pdf23,pdf31,pdf32,pdf33,pdf,lnf, Kt11,Kt12,Kt13,Kt21,Kt22,Kt23,Kt31,Kt32,Kt33, beta_tt11,beta_tt12,beta_tt13,beta_tt21,beta_tt22,beta_tt23,beta_tt31,beta_tt32,beta_tt33, P_tt11,P_tt12,P_tt13,P_tt21,P_tt22,P_tt23,P_tt31,P_tt32,P_tt33, beta_tt1, beta_tt2, beta_tt3,P_tt1,P_tt2,P_tt3,bb2,bb1,bb3,pprl,PR_TR,PPR_tl,J_pdf, cycle,trend,ss2mat,ss3mat,SD_resid,eta_mat,f_mat,f_tl,eta_tl,beta_tt ,cor1,cor2,cor3,cov1,cov2,cov3,mu1,mu2,mu3; prmtr=trans(prmtr1); aa1=prmtr[1]; aa2=prmtr[2]; aa3=prmtr[3]; bb1=prmtr[4]; bb2=prmtr[5]; bb3=prmtr[6]; ssig2v1=prmtr[7]^2; ssig2v2=prmtr[8]^2; ssig2v3=prmtr[9]^2; ssig2w1=prmtr[10]^2; ssig2w2=prmtr[11]^2; ssig2w3=prmtr[12]^2; p11=prmtr[13];@Prb[ S(t)=1 | S(t-1)=1]@ p22=prmtr[14];@Prb[ S(t)=0 | S(t-1)=0]@ cycle=zeros(t,1); trend=zeros(t,1); ss2mat=zeros(t,2); ss3mat=zeros(t,2); SD_resid=zeros(t,1); eta_mat=zeros(t,1); f_mat=zeros(t,1); cor1=prmtr[15]; cor2=prmtr[16]; cor3=prmtr[17]; mu1=prmtr[18]; mu2=prmtr[19]; mu3=prmtr[20]; @ cor1=1; cor2=1; cor3=1; @ cov1=cor1*sqrt(ssig2v1*ssig2w1); cov2=cor2*sqrt(ssig2v2*ssig2w2); cov3=cor3*sqrt(ssig2v3*ssig2w3); @ F_j @ F1=(1~ 0 ~ 0 )| (0~aa1~bb1)| (0~ 1 ~ 0 ); F2=(1~ 0 ~ 0 )| (0~aa2~bb2)| (0~ 1 ~ 0 ); F3=(1~ 0 ~ 0 )| (0~aa3~bb3)| (0~ 1 ~ 0 ); @ Q_j @ Q1=( ssig2v1 ~ cov1 ~ 0 )| ( cov1 ~ ssig2w1 ~ 0 )| ( 0 ~ 0 ~ 0 ); Q2=( ssig2v2 ~ cov2 ~ 0 )| ( cov2 ~ ssig2w2 ~ 0 )| ( 0 ~ 0 ~ 0 ); Q3=( ssig2v3 ~ cov3 ~ 0 )| ( cov3 ~ ssig2w3 ~ 0 )| ( 0 ~ 0 ~ 0 ); R1=zeros(n,n); R2=zeros(n,n); R3=zeros(n,n); @ n by n @ beta_ll1=ymat[1]|0|0; beta_ll2=ymat[1]|0|0; beta_ll3=ymat[1]|0|0; p_ll1=500*eye(k); p_ll2=500*eye(k); p_ll3=500*eye(k);@ p_ll1=5*eye(k); p_ll2=5*eye(k); p_ll3=5*eye(k);@ p12=1-p11; p13=0; p21=0; p23=1-p22; p31=0; p32=0; p33=1; PR_TR= p11 ~ p12 ~ p13| p21 ~ p22 ~ p23| p31 ~ p32 ~ p33; PR_TR=vec(PR_TR); @ 9 by 1 @ @ At t=0, to calculate conditional probability St=1 , St=2 and St=3 @ pprl1_=0.99998; @ Unconditional Probability for St=1 @ pprl2_=0.00001; @ Unconditional Probability for St=2 @ pprl3_=0.00001; @ Unconditional Probability for St=3 @ pprl=pprl1_|pprl2_|pprl3_; @ 3 by 1 @ pprl=vec(pprl~pprl~pprl); @ 9 by 1 @ lnL=0; itr=1; do until itr>t; @ == Step 1 : Kalman Filter == @ @ beta_tl(i,j), S(t-1)=i, S(t)=j , k by 1 @ beta_tl11=mu1+ F1*beta_ll1; @ S(t-1)=1, S(t)=1 @ beta_tl21=mu1+ F1*beta_ll2; @ S(t-1)=2, S(t)=1 @ beta_tl31=mu1+ F1*beta_ll3; @ S(t-1)=3, S(t)=1 @ beta_tl12=mu2+ F2*beta_ll1; @ S(t-1)=1, S(t)=2 @ beta_tl22= mu2+F2*beta_ll2; @ S(t-1)=2, S(t)=2 @ beta_tl32=mu2+ F2*beta_ll3; @ S(t-1)=3, S(t)=2 @ beta_tl13= mu3+F3*beta_ll1; @ S(t-1)=1, S(t)=3 @ beta_tl23=mu3+ F3*beta_ll2; @ S(t-1)=2, S(t)=3 @ beta_tl33=mu3+ F3*beta_ll3; @ S(t-1)=3, S(t)=3 @ @ P_tl(i,j), S(t-1)=i, S(t)=j , k by k @ P_tl11=F1*P_ll1*F1' +Q1; P_tl21=F1*P_ll2*F1' +Q1; P_tl31=F1*P_ll3*F1' +Q1; P_tl12=F2*P_ll1*F2' +Q2; P_tl22=F2*P_ll2*F2' +Q2; P_tl32=F2*P_ll3*F2' +Q2; P_tl13=F3*P_ll1*F3' +Q3; P_tl23=F3*P_ll2*F3' +Q3; P_tl33=F3*P_ll3*F3' +Q3; @ H_t( j ), S(t)=j @ H_t1=1~1~0; @ 1 by k @ H_t2=1~1~0; H_t3=1~1~0; @ eta_tl(i,j), S(t-1)=i, S(t)=j , n by n @ eta_tl11=ymat[itr]-H_t1*beta_tl11 ; eta_tl21=ymat[itr]-H_t1*beta_tl21 ; eta_tl31=ymat[itr]-H_t1*beta_tl31 ; eta_tl12=ymat[itr]-H_t2*beta_tl12 ; eta_tl22=ymat[itr]-H_t2*beta_tl22 ; eta_tl32=ymat[itr]-H_t2*beta_tl32 ; eta_tl13=ymat[itr]-H_t3*beta_tl13 ; eta_tl23=ymat[itr]-H_t3*beta_tl23 ; eta_tl33=ymat[itr]-H_t3*beta_tl33 ; @ f_tl(i,j), S(t-1)=i, S(t)=j , n by n @ f_tl11=H_t1*P_tl11*H_t1'- R1; f_tl21=H_t1*P_tl21*H_t1'- R1; f_tl31=H_t1*P_tl31*H_t1'- R1; f_tl12=H_t2*P_tl12*H_t2'- R2; f_tl22=H_t2*P_tl22*H_t2'- R2; f_tl32=H_t2*P_tl32*H_t2'- R2; f_tl13=H_t3*P_tl13*H_t3'- R3; f_tl23=H_t3*P_tl23*H_t3'- R3; f_tl33=H_t3*P_tl33*H_t3'- R3; @ Kt(i,j), S(t-1)=i, S(t)=j , k by n @ Kt11=P_tl11*H_t1'*inv(f_tl11); Kt21=P_tl21*H_t1'*inv(f_tl21); Kt31=P_tl31*H_t1'*inv(f_tl31); Kt12=P_tl12*H_t2'*inv(f_tl12); Kt22=P_tl22*H_t2'*inv(f_tl22); Kt32=P_tl32*H_t2'*inv(f_tl32); Kt13=P_tl13*H_t3'*inv(f_tl13); Kt23=P_tl23*H_t3'*inv(f_tl23); Kt33=P_tl33*H_t3'*inv(f_tl33); @ == Updating == @ @ beta_tl(i,j), S(t-1)=i, S(t)=j , k by 1 @ beta_tt11= beta_tl11 + Kt11*eta_tl11; beta_tt21= beta_tl21 + Kt21*eta_tl21; beta_tt31= beta_tl31 + Kt31*eta_tl31; beta_tt12= beta_tl12 + Kt12*eta_tl12; beta_tt22= beta_tl22 + Kt22*eta_tl22; beta_tt32= beta_tl32 + Kt32*eta_tl32; beta_tt13= beta_tl13 + Kt13*eta_tl13; beta_tt23= beta_tl23 + Kt23*eta_tl23; beta_tt33= beta_tl33 + Kt33*eta_tl33; @ P_tl(i,j), S(t-1)=i, S(t)=j , k by k @ P_tt11=(eye(k)-P_tl11*H_t1'*inv(f_tl11)*H_t1)*P_tl11; P_tt21=(eye(k)-P_tl21*H_t1'*inv(f_tl21)*H_t1)*P_tl21; P_tt31=(eye(k)-P_tl31*H_t1'*inv(f_tl31)*H_t1)*P_tl31; P_tt12=(eye(k)-P_tl12*H_t2'*inv(f_tl12)*H_t2)*P_tl12; P_tt22=(eye(k)-P_tl22*H_t2'*inv(f_tl22)*H_t2)*P_tl22; P_tt32=(eye(k)-P_tl32*H_t2'*inv(f_tl32)*H_t2)*P_tl32; P_tt13=(eye(k)-P_tl13*H_t3'*inv(f_tl13)*H_t3)*P_tl13; P_tt23=(eye(k)-P_tl23*H_t3'*inv(f_tl23)*H_t3)*P_tl23; P_tt33=(eye(k)-P_tl33*H_t3'*inv(f_tl33)*H_t3)*P_tl33; @ == Step 2 : Hamilton Filter == @ PPR_tl=pprl.* PR_TR; @ Pr[S(t),S(t-1) l I(t-1)], 9 by 1 @ ss2mat[itr,1]=sumc(PPR_tl[4:6]);@ Pr[S(t) l (t-1)], 9 by 1 @ ss3mat[itr,1]=sumc(PPR_tl[7:9]); pdf11=(1/sqrt(2*pi*det(f_tl11)))*exp(-0.5*eta_tl11'inv(f_tl11)*eta_tl11); pdf21=(1/sqrt(2*pi*det(f_tl21)))*exp(-0.5*eta_tl21'inv(f_tl21)*eta_tl21); pdf31=(1/sqrt(2*pi*det(f_tl31)))*exp(-0.5*eta_tl31'inv(f_tl31)*eta_tl31); pdf12=(1/sqrt(2*pi*det(f_tl12)))*exp(-0.5*eta_tl12'inv(f_tl12)*eta_tl12); pdf22=(1/sqrt(2*pi*det(f_tl22)))*exp(-0.5*eta_tl22'inv(f_tl22)*eta_tl22); pdf32=(1/sqrt(2*pi*det(f_tl32)))*exp(-0.5*eta_tl32'inv(f_tl32)*eta_tl32); pdf13=(1/sqrt(2*pi*det(f_tl13)))*exp(-0.5*eta_tl13'inv(f_tl13)*eta_tl13); pdf23=(1/sqrt(2*pi*det(f_tl23)))*exp(-0.5*eta_tl23'inv(f_tl23)*eta_tl23); pdf33=(1/sqrt(2*pi*det(f_tl33)))*exp(-0.5*eta_tl33'inv(f_tl33)*eta_tl33); pdf=pdf11|pdf21|pdf31|pdf12|pdf22|pdf32|pdf13|pdf23|pdf33; pdf=pdf.*PPR_tl; eta_tl=eta_tl11|eta_tl21|eta_tl31|eta_tl12|eta_tl22|eta_tl32|eta_tl13|eta_tl23|eta_tl33; eta_mat[itr]=sumc(eta_tl.*PPR_tl); @ Prediction error @ f_tl=f_tl11|f_tl21|f_tl31|f_tl12|f_tl22|f_tl32|f_tl13|f_tl23|f_tl33; f_mat[itr]=sumc(f_tl.*PPR_tl); @ Variance of prediction error @ SD_resid[itr]=eta_mat[itr]/sqrt(f_mat[itr]); @ SD_resid[itr]=sumc(eta_tl.*ppr_tl./sqrt(f_tl));@ J_pdf=sumc(pdf); if itr>start; lnf=ln(J_pdf); else; lnf=0; endif; lnL=lnL+lnf; pprl1_=sumc(pdf[1:3])/J_pdf; pprl2_=sumc(pdf[4:6])/J_pdf; pprl3_=sumc(pdf[7:9])/J_pdf; @ == Step 3 : Approximations == @ @ beta_tt(j), S(t) = j, k by 1 @ beta_tt1=(beta_tt11*(pdf[1]/j_pdf)+beta_tt21*(pdf[2]/j_pdf)+beta_tt31*(pdf[3]/j_pdf))/pprl1_; @ beta_tt(1) @ beta_tt2=(beta_tt12*(pdf[4]/j_pdf)+beta_tt22*(pdf[5]/j_pdf)+beta_tt32*(pdf[6]/j_pdf))/pprl2_; @ beta_tt(1) @ beta_tt3=(beta_tt13*(pdf[7]/j_pdf)+beta_tt23*(pdf[8]/j_pdf)+beta_tt33*(pdf[9]/j_pdf))/pprl3_; @ beta_tt(1) @ @ P_tt(j), S(t) = j, k by k @ P_tt1=( (pdf[1]/J_pdf)*(P_tt11+(beta_tt1-beta_tt11)*(beta_tt1-beta_tt11)') +(pdf[2]/J_pdf)*(P_tt21+(beta_tt1-beta_tt21)*(beta_tt1-beta_tt21)') +(pdf[3]/J_pdf)*(P_tt31+(beta_tt1-beta_tt31)*(beta_tt1-beta_tt31)') )/pprl1_; @ P_tt(1) @ P_tt2=( (pdf[4]/J_pdf)*(P_tt12+(beta_tt2-beta_tt12)*(beta_tt2-beta_tt12)') +(pdf[5]/J_pdf)*(P_tt22+(beta_tt2-beta_tt22)*(beta_tt2-beta_tt22)') +(pdf[6]/J_pdf)*(P_tt32+(beta_tt2-beta_tt32)*(beta_tt2-beta_tt32)') )/pprl2_; @ P_tt(1) @ P_tt3=( (pdf[7]/J_pdf)*(P_tt13+(beta_tt3-beta_tt13)*(beta_tt3-beta_tt13)') +(pdf[8]/J_pdf)*(P_tt23+(beta_tt3-beta_tt23)*(beta_tt3-beta_tt23)') +(pdf[9]/J_pdf)*(P_tt33+(beta_tt3-beta_tt33)*(beta_tt3-beta_tt33)') )/pprl3_; @ P_tt(1) @ pprl=pprl1_|pprl2_|pprl3_; @ 3 by 1 @ pprl=vec(pprl~pprl~pprl); @ 9 by 1 @ beta_ll1=beta_tt1; beta_ll2= beta_tt2; beta_ll3= beta_tt3; beta_tt= pprl1_*beta_tt1 + pprl2_*beta_tt2+pprl3_*beta_tt3; trend[itr]=beta_tt[1]; cycle[itr]=beta_tt[2]; ss2mat[itr,2]=pprl2_; ss3mat[itr,2]=pprl3_; P_ll1= P_tt1;P_ll2= P_tt2;P_ll3= P_tt3; itr=itr+1; endo; retp(trend[1:t], cycle[1:t], ss2mat[1:t,.],ss3mat[1:t,.],SD_resid[1:t,.],eta_mat[1:t,.],f_mat[1:t,.]); endp; /*===================== Procedure 3: Kim's Smoothing =====================*/ proc kim(prb_h,prb_tl); local pr_sm1,pr_sm2,pr_sm3, j_iter,pr_sm33,pr_sm12,pr_sm21,pr_sm22,pr_sm11,pr_sm23, pr_tt1,pr_tt2,pr_tt3,pr_tl1,pr_tl2,pr_tl3,i; pr_tt1=prb_h[.,1]; pr_tt2=prb_h[.,2]; pr_tt3=prb_h[.,3]; pr_tl1=prb_tl[.,1]; pr_tl2=prb_tl[.,2]; pr_tl3=prb_tl[.,3]; pr_sm1=pr_tt1; @ pr_sm0 will contain Pr[S_t|Y_T]@ pr_sm2=pr_tt2; pr_sm3=pr_tt3; j_iter=T-1; do until j_iter < 1; @The followings are P[S_t, S_t+1|Y_T] @ pr_sm11=pr_sm1[j_iter+1,1]*p11* pr_tt1[j_iter,1]/ pr_tl1[j_iter+1,1]; pr_sm12=pr_sm2[j_iter+1,1]*p12* pr_tt1[j_iter,1]/ pr_tl2[j_iter+1,1]; pr_sm22=pr_sm2[j_iter+1,1]*p22* pr_tt2[j_iter,1]/ pr_tl2[j_iter+1,1]; pr_sm23=pr_sm3[j_iter+1,1]*p23* pr_tt2[j_iter,1]/ pr_tl3[j_iter+1,1]; pr_sm33=pr_sm3[j_iter+1,1]*p33* pr_tt3[j_iter,1]/ pr_tl3[j_iter+1,1]; pr_sm1[j_iter,1]=pr_sm11+pr_sm12; pr_sm2[j_iter,1]=pr_sm22+pr_sm23; pr_sm3[j_iter,1]=pr_sm33; j_iter=j_iter -1; endo; RETP(pr_sm2~pr_sm3); @This proc returns Pr[S_t=1 & 2 |Y_T]@ endp; /*=================== Transformation ===================*/ proc trans(prmtr1); local prmtr,lam11,lam12,lam21,lam22,lam31,lam32; prmtr=prmtr1; /* lam11=(prmtr1[1])/(abs(prmtr1[1])+1); lam12=(prmtr1[4])/(abs(prmtr1[4])+1); lam21=(prmtr1[2])/(abs(prmtr1[2])+1); lam22=(prmtr1[5])/(abs(prmtr1[5])+1); lam31=(prmtr1[3])/(abs(prmtr1[3])+1); lam32=(prmtr1[6])/(abs(prmtr1[6])+1); prmtr[1]=lam11+lam12; prmtr[2]=lam21+lam22; prmtr[3]=lam31+lam32; prmtr[4]=-lam11*lam12; prmtr[5]=-lam21*lam22; prmtr[6]=-lam31*lam32; */ prmtr[7:12]=sqrt(exp(prmtr1[7:12]));@sig2@ prmtr[13:14]=exp(prmtr1[13:14])./(exp(prmtr1[13:14])+1); @ p11,p22 @ prmtr[15:17]=100*(prmtr[15:17])./(abs(prmtr[15:17])*100+1); retp(prmtr); endp; proc trans0(prmtr1); local prmtr,lam11,lam12,lam21,lam22,lam31,lam32; prmtr=prmtr1; /* lam11=(prmtr1[1])/(abs(prmtr1[1])+1); lam12=(prmtr1[4])/(abs(prmtr1[4])+1); lam21=(prmtr1[2])/(abs(prmtr1[2])+1); lam22=(prmtr1[5])/(abs(prmtr1[5])+1); lam31=(prmtr1[3])/(abs(prmtr1[3])+1); lam32=(prmtr1[6])/(abs(prmtr1[6])+1); prmtr[1]=lam11+lam12; prmtr[2]=lam21+lam22; prmtr[3]=lam31+lam32; prmtr[4]=-lam11*lam12; prmtr[5]=-lam21*lam22; prmtr[6]=-lam31*lam32; */ prmtr[7:10]=sqrt(exp(prmtr1[7:10]));@sig2@ prmtr[11:12]=exp(prmtr1[11:12])./(exp(prmtr1[11:12])+1); @ p11,p22 @ prmtr[13:15]=100*(prmtr[13:15])./(abs(prmtr[13:15])*100+1); retp(prmtr); endp; @ Trans for Model Transformation @ proc arma(prmtr1); local prmtr,a1,a2,a3,b1,b2,b3,sig2v1,sig2v2,sig2v3,sig2w1,sig2w2,sig2w3, phi11,phi21,phi12,phi22,phi13,phi23,sig2e1,sig2e2,sig2e3, theta11,theta21,theta12,theta22,theta13,theta23,cov1,cov2,cov3,cor1,cor2,cor3,A,B,C,coef1,rt1,coef2,coef3,rt2,rt3; @prmtr1=trans(prmtr1);@ prmtr=zeros(15,1); a1=prmtr1[1]; a2=prmtr1[2]; a3=prmtr1[3]; b1=prmtr1[4]; b2=prmtr1[5]; b3=prmtr1[6]; sig2v1=prmtr1[7]^2; sig2v2=prmtr1[8]^2; sig2v3=prmtr1[9]^2; sig2w1=prmtr1[10]^2; sig2w2=prmtr1[11]^2; sig2w3=prmtr1[12]^2; cor1=prmtr1[15]; cor2=prmtr1[16]; cor3=prmtr1[17]; cov1=cor1*sqrt(sig2v1*sig2w1); cov2=cor2*sqrt(sig2v2*sig2w2); cov3=cor3*sqrt(sig2v3*sig2w3); phi11=a1; phi21=b1; A=sig2v1*(1+phi11^2+phi21^2)+2*sig2w1+(1+phi11)*cov1; B=phi11*(phi21-1)*sig2v1-sig2w1+cov1*(phi21-phi11-1); C=-phi21*(sig2v1+cov1); coef1=1|(2*C-A)|(2*C^2+B^2-2*A*C)|(2*C^3-A*C^2)|(C^4); rt1=polyroot(coef1); @"rt1"; rt1;@ sig2e1=abs(rt1[1]); theta11=B/(sig2e1+C); theta21=C/sig2e1; phi12=a2; phi22=b2; A=sig2v2*(1+phi12^2+phi22^2)+2*sig2w2+(1+phi12)*cov2; B=phi12*(phi22-1)*sig2v2-sig2w2+cov2*(phi22-phi12-1); C=-phi22*(sig2v2+cov2); coef2=1|(2*C-A)|(2*C^2+B^2-2*A*C)|(2*C^3-A*C^2)|(C^4); rt2=polyroot(coef2); @"rt2"; rt2;@ sig2e2=abs(rt2[1]); theta12=B/(sig2e2+C); theta22=C/sig2e2; phi13=a3; phi23=b3; A=sig2v3*(1+phi13^2+phi23^2)+2*sig2w3+(1+phi13)*cov3; B=phi13*(phi23-1)*sig2v3-sig2w3+cov3*(phi23-phi13-1); C=-phi23*(sig2v3+cov3); coef3=1|(2*C-A)|(2*C^2+B^2-2*A*C)|(2*C^3-A*C^2)|(C^4); rt3=polyroot(coef3); @"rt3"; rt3;@ sig2e3=abs(rt3[1]); theta13=B/(sig2e3+C); theta23=C/sig2e3; prmtr[1]=phi11; prmtr[2]=phi21; prmtr[3]=theta11; prmtr[4]=theta21; prmtr[5]=sig2e1; prmtr[6]=phi12; prmtr[7]=phi22; prmtr[8]=theta12; prmtr[9]=theta22; prmtr[10]=sig2e2; prmtr[11]=phi13; prmtr[12]=phi23; prmtr[13]=theta13; prmtr[14]=theta23; prmtr[15]=sig2e3; retp(prmtr); endp; proc duration(prmtr); local out; out=1/(1-prmtr); retp(out); endp;