/* IDENTIFYING BUSINESS CYCLE PHASES WITH COINCIDENT ECONOMIC VARIABLES: A GIBBS SAMPLING APPROACH TO DYNAMIC FACTOR MODEL WITH MARKOV-SWITCHING PROGRAM WRITTEN BY CHANG-JIN KIM DEPT. OF ECONOMICS KOREA UNIVERSITY, SEOUL 136-701 */ new; #lineson; library pgraph; clear counter; clear totcount; format /m1 /rd 14,8; LOAD DATA[517,6]=FULLDTA.PRN; DATA1=DATA[.,2]; DATA2=DATA[.,3]; DATA3=DATA[.,4]; DATA4=DATA[.,5]; lndata1=ln(data1)*100; lndata2=ln(data2)*100; lndata3=ln(data3)*100; lndata4=ln(data4)*100; t0=rows(lndata1); y1= lndata1[2:t0,1]-lndata1[1:t0-1,1]; y2= lndata2[2:t0,1]-lndata2[1:t0-1,1]; y3= lndata3[2:t0,1]-lndata3[1:t0-1,1]; y4= lndata4[2:t0,1]-lndata4[1:t0-1,1]; @ yy0=y1/stdc(y1)~y2/stdc(y2)~y3/stdc(y3)~y4/stdc(y4);@ yy0=y1~y2~y3~y4; begin_t=85; @85---59.02@ @96---60.01@ end_t=516; @1952.2--1995.01@ yy1=yy0[BEGIN_T:END_T,.]; @4xT@ T=rows(yy1); seed=198011; rndns(1,1,seed); @ yyy= yy1'./ (stdc(yy1).*ones(4,T));@ yy_m=meanc(yy1); yy= yy1' - yy_m; t=cols(yy); @================= PRIOR PARAMETERS ====================@ U1_01_=1; U1_00_=1; @FOR A1TT@ U1_10_=1; U1_11_=1; @FOR B1TT@ R0_=EYE(2)/4; T0_=0|0; @FOR PHI1,PHI2 OF COMMON COMPONENT@ R0_M=EYE(2)/4; @For MU0, MU1 OF COMMON COMPONENT@ T0_M=0|0; R0_V=EYE(2)/4; T0_V=0|0; @FOR PSI1,PSI2 OF INDIVIDUAL COMPONENT@ D0_=0; V0_=0; @FOR SIG2_V OF INDIVIDUAL COMPONENT@ R00_=1/4; T00_=0; @FOR LAMBDA_i OF INDIVIDUAL COMPONENT@ R00_4=EYE(4)/1; T00_4=0|0|0|0; @FOR LAMDA'S OF EMPLOYMENT EQN@ @========INITIAL VALUES FOR THE GIBBS SAMPLER ===========@ A1TT=1-0.9; B1TT=1-0.97; PHI1TT=0; PHI2TT=0; MU0TT=-2.0; MU1TT=2.5; LAMD1TT=0.5; PSI11TT=0; PSI12TT=0; SIG21TT=0.2; LAMD2TT=0.5; PSI21TT=0; PSI22TT=0; SIG22TT=0.2; LAMD3TT=0.5; PSI31TT=0; PSI32TT=0; SIG23TT=0.2; LAMD4TT=0.5; PSI41TT=0; PSI42TT=0; SIG24TT=0.2; LAM41TT=0; LAM42TT=0;LAM43TT=0; STT=markov(a1tt,b1tt,1,T); stt[1:2,1]=ONES(2,1); @Check this@ @===== DESIGN FOR THE GIBBS SAMPLER ===========@ N0 = 2000; @NUMBER OF DRAWS TO LEAVE OUT@ MM = 8000; @NUMBER OF DRAWS@ LL = 5; @EVERY L-TH VALUE IS USED@ CAPN = N0 + MM; @TOTAL NUMBER OF DRAWS@ @===== STORAGE SPACES===========================@ PSI11MM = {}; PSI12MM = {}; PSI21MM = {}; PSI22MM = {}; PSI31MM = {}; PSI32MM = {}; PSI41MM = {}; PSI42MM = {}; SIG21MM = {}; SIG22MM = {}; SIG23MM = {}; SIG24MM = {}; LAMD1MM = {}; LAMD2MM = {}; LAMD3MM = {}; LAMD4MM = {}; LAM41MM ={}; LAM42MM={}; LAM43MM={}; PHI1MM = {}; PHI2MM = {}; MU0MM = {}; MU1MM = {}; A1MM = {}; B1MM = {}; MUMM ={}; STTMM=ZEROS(T,1); ZTTMM=ZEROS(T,1); CCIMM=ZEROS(T,1); dltmm={}; STTMM2=ZEROS(T,1); CCIMM2=ZEROS(T,1); swt_mat=zeros(t,4); ITR=1; DO WHILE ITR LE CAPN; ITR;;"th iteration......."; @=================START SAMPLING================@ ZTT = GEN_ZT; STT = GEN_ST; swtmp=sw_cnt(stt); {PSI11TT,PSI12TT,SIG21TT}=GEN_PSI(1,LAMD1TT,SIG21TT); {PSI21TT,PSI22TT,SIG22TT}=GEN_PSI(2,LAMD2TT,SIG22TT); {PSI31TT,PSI32TT,SIG23TT}=GEN_PSI(3,LAMD3TT,SIG23TT); {PSI41TT,PSI42TT,SIG24TT}=GEN_PSI4(4,LAMD4TT,SIG24TT); LAMD1TT = GEN_LMDA(1,PSI11TT,PSI12TT,SIG21TT); LAMD2TT = GEN_LMDA(2,PSI21TT,PSI22TT,SIG22TT); LAMD3TT = GEN_LMDA(3,PSI31TT,PSI32TT,SIG23TT); {LAMD4TT,LAM41TT,LAM42TT,LAM43TT} = GEN_LMD4(4,PSI41TT,PSI42TT,SIG24TT); {phi1tt,phi2tt,mu0tt,mu1tt} = GEN_PHI; locate 10,1; tranmat = switchg(stt[5:T,1]+1,1|2); print "this is tranmat" tranmat; "mu";; mu0tt;;mu1tt; "lmdt";;lamd1tt;;lamd2tt; " ";;lamd3tt;;lamd4tt; A1TT = rndb(tranmat[1,2]+u1_01_,tranmat[1,1]+u1_00_); B1TT = rndb(tranmat[2,1]+u1_10_,tranmat[2,2]+u1_11_); {CCI,dlt} = GEN_CCI; @xy(seqa(1,1,rows(cci[3:t])),cci[3:t]);@ @==========END OF ONE ITERATION ========================@ PSI11MM = PSI11MM~PSI11TT; PSI12MM = PSI12MM~PSI12TT; PSI21MM = PSI21MM~PSI21TT; PSI22MM = PSI22MM~PSI22TT; PSI31MM = PSI31MM~PSI31TT; PSI32MM = PSI32MM~PSI32TT; PSI41MM = PSI41MM~PSI41TT; PSI42MM = PSI42MM~PSI42TT; SIG21MM = SIG21MM~SIG21TT; SIG22MM = SIG22MM~SIG22TT; SIG23MM = SIG23MM~SIG23TT; SIG24MM = SIG24MM~SIG24TT; LAMD1MM = LAMD1MM~LAMD1TT; LAMD2MM = LAMD2MM~LAMD2TT; LAMD3MM = LAMD3MM~LAMD3TT; LAMD4MM = LAMD4MM~LAMD4TT; LAM41MM=LAM41MM~LAM41TT; LAM42MM=LAM42MM~LAM42TT; LAM43MM=LAM43MM~LAM43TT; PHI1MM = PHI1MM~PHI1TT; PHI2MM = PHI2MM~PHI2TT; MU0MM = MU0MM~MU0TT; MU1MM = MU1MM~MU1TT; MUMM=MUMM~(MU0TT+MU1TT); A1MM = A1MM~(1-A1TT) ; B1MM = B1MM~(1-B1TT); dltmm=dltmm~dlt; if itr>n0; STTMM=STTMM+STT; ZTTMM=ZTTMM+ZTT; CCIMM=CCIMM+CCI; STTMM2=STTMM2+STT.*STT; CCIMM2=CCIMM2+ CCI.*CCI; swt_mat=swt_mat+swtmp; @00,01,10,11@ ENDIF; @if itr>100; xy(seqa(1,1,rows(cci)),cci); endif; @ ITR=ITR+1; ENDO; STTMM=STTMM/MM; ZTTMM=ZTTMM/MM; CCIMM=CCIMM/MM; FNL_MAT=A1MM'~B1MM'~PHI1MM'~PHI2MM'~MU0MM'~MU1MM'~ LAMD1MM'~PSI11MM'~PSI12MM'~SIG21MM'~ LAMD2MM'~PSI21MM'~PSI22MM'~SIG22MM'~ LAMD3MM'~PSI31MM'~PSI32MM'~SIG23MM'~ LAMD4MM'~PSI41MM'~PSI42MM'~SIG24MM'~dltmm'~LAM41MM'~ LAM42MM'~LAM43MM'~MUMM'; SAVE FNL_MAT,STTMM,ZTTMM,CCIMM,STTMM2,CCIMM2,swt_mat; @======== calculation of posterior expectations and std deviations@ fnl_mat=fnl_mat[n0+1:rows(fnl_mat),.]; INDX = ZEROS(MM,1); INDX[1] = 1; J = 1; DO WHILE J LE MM; INDX[J] = 1; J = J + LL; ENDO; SORT_OUT=ZEROS(cols(fnl_mat),3); I_NDX=1; DO UNTIL I_NDX>COLS(FNL_MAT); tmpm1 = selif(FNL_MAT[.,I_NDX],indx); MN_OUT = meanc(tmpm1); STD_OUT = stdc(tmpm1); MED_OUT = median(tmpm1); u2out = MN_OUT~STD_OUT~MED_OUT; SORT_OUT[I_NDX,.]=U2OUT; I_NDX=I_NDX+1; ENDO; OUTPUT FILE=MMAIN_4.OUT RESET; "FROM 1st TO 22nd ROWS.............."; " FNL_MAT=A1MM'~B1MM'~PHI1MM'~PHI2MM'~MU0MM'~MU1MM'~ LAMD1MM'~PSI11MM'~PSI12MM'~SIG21MM'~ LAMD2MM'~PSI21MM'~PSI22MM'~SIG22MM'~ LAMD3MM'~PSI31MM'~PSI32MM'~SIG23MM'~ LAMD4MM'~PSI41MM'~PSI42MM'~SIG24MM'~dltmm' ~LAM41MM~LAM42MM~LAM43MM; "; SORT_OUT; OUTPUT OFF; END; @========================================================@ @ A PROCEDURE THAT GENERATES PSI_1,PSI_2,SIG2_V @ @========================================================@ PROC(3) = GEN_PSI(NDX,LAMDA,SIG2_V); local ystarr,ystar,xstar,v,psi,c,accept,psi_g,root,rootmod,nn,d,t2, sig2_v_g,coef,tstar; YSTARR=YY[NDX,1:T]'-LAMDA*ZTT[1:T,1]; YSTAR=YSTARR[5:T,1]; XSTAR=YSTARR[4:T-1,1]~YSTARR[3:T-2,1]; tstar=rows(ystar); V = invpd(R0_v + SIG2_V^(-1)*XSTAR'XSTAR); PSI = V*(R0_v*t0_v + SIG2_V^(-1)*XSTAR'YSTAR); C = chol(V); ACCEPT = 0; DO WHILE ACCEPT ==0; PSI_G = PSI + C'rndn(2,1); COEF = -REV(PSI_G)|1; ROOT = POLYROOT(COEF); ROOTMOD = ABS(ROOT); IF MINC(ROOTMOD) GE 1.0001; ACCEPT = 1; ELSE; COUNTER = COUNTER + 1; ACCEPT = 0; ENDIF; ENDO; nn = tstar + v0_; d = d0_ + (YSTAR-XSTAR*PSI_G)'(YSTAR-XSTAR*PSI_G); c = rndc(nn); t2 = c/d; SIG2_V_G = 1/t2; RETP(PSI_G[1,1],PSI_G[2,1],SIG2_V_G); ENDP; PROC(3) = GEN_PSI4(NDX,LAMDA,SIG2_V); local ystarr,ystar,xstar,v,psi,c,accept,psi_g,root,rootmod,nn,d,t2, sig2_v_g,coef,tstar; YSTARR=YY[NDX,4:T]'-LAMDA*ZTT[4:T,1]- LAM41TT*ZTT[3:T-1,1] - LAM42TT*ZTT[2:T-2,1]-LAM43TT*ZTT[1:T-3,1]; YSTAR=YSTARR[5:T-3,1]; XSTAR=YSTARR[4:T-3-1,1]~YSTARR[3:T-3-2,1]; tstar=rows(ystar); V = invpd(R0_v + SIG2_V^(-1)*XSTAR'XSTAR); PSI = V*(R0_v*t0_v + SIG2_V^(-1)*XSTAR'YSTAR); C = chol(V); ACCEPT = 0; DO WHILE ACCEPT ==0; PSI_G = PSI + C'rndn(2,1); COEF = -REV(PSI_G)|1; ROOT = POLYROOT(COEF); ROOTMOD = ABS(ROOT); IF MINC(ROOTMOD) GE 1.0001; ACCEPT = 1; ELSE; COUNTER = COUNTER + 1; ACCEPT = 0; ENDIF; ENDO; nn = tstar + v0_; d = d0_ + (YSTAR-XSTAR*PSI_G)'(YSTAR-XSTAR*PSI_G); c = rndc(nn); t2 = c/d; SIG2_V_G = 1/t2; RETP(PSI_G[1,1],PSI_G[2,1],SIG2_V_G); ENDP; @========================================================@ @ A PROCEDURE THAT GENERATES LAMBDA @ @========================================================@ PROC GEN_LMDA(NDX,PSI_1,PSI_2,SIG2_V); local ystar,xstar,v,lamda,c,lamda_g; YSTAR=YY[NDX,5:T]'-PSI_1*YY[NDX,4:T-1]'-PSI_2*YY[NDX,3:T-2]'; XSTAR=ZTT[5:T,1] -PSI_1*ZTT[4:T-1,1] -PSI_2*ZTT[3:T-2,1]; V = invpd(R00_ + SIG2_V^(-1)*XSTAR'XSTAR); LAMDA = V*(R00_*t00_ + SIG2_V^(-1)*XSTAR'YSTAR); C = chol(V); LAMDA_G = LAMDA + C'rndn(1,1); RETP(LAMDA_G); ENDP; PROC(4)= GEN_LMD4(NDX,PSI_1,PSI_2,SIG2_V); local ystar,xstar,v,lamda,c,lamda_g; YSTAR=YY[NDX,8:T]'-PSI_1*YY[NDX,7:T-1]'-PSI_2*YY[NDX,6:T-2]'; XSTAR=(ZTT[8:T,1] -PSI_1*ZTT[7:T-1,1] -PSI_2*ZTT[6:T-2,1])~ (ZTT[7:T-1,1] -PSI_1*ZTT[6:T-2,1] -PSI_2*ZTT[5:T-3,1])~ (ZTT[6:T-2,1] -PSI_1*ZTT[5:T-3,1] -PSI_2*ZTT[4:T-4,1])~ (ZTT[5:T-3,1] -PSI_1*ZTT[4:T-4,1] -PSI_2*ZTT[3:T-5,1]); V = invpd(R00_4 + SIG2_V^(-1)*XSTAR'XSTAR); LAMDA = V*(R00_4*t00_4 + SIG2_V^(-1)*XSTAR'YSTAR); C = chol(V); LAMDA_G = LAMDA + C'rndn(4,1); RETP(LAMDA_G[1,1],LAMDA_G[2,1],LAMDA_G[3,1],LAMDA_G[4,1]); ENDP; @========================================================@ @ A PROCEDURE THAT GENERATES mu0,mu1 and phi1,phi2 @ @========================================================@ PROC (4) = GEN_PHI; local ystar,xstar,v,c,accept,coef,root,rootmod, tstar,YSTARR,PHI_G,MU_G,PHI,MU; YSTARR=ZTT-MU0TT*ONES(T,1)-MU1TT*STT; YSTAR=YSTARR[5:T,1]; XSTAR=YSTARR[4:T-1,1]~YSTARR[3:T-2,1]; V = invpd(R0_ + XSTAR'XSTAR); PHI = V*(R0_*t0_ + XSTAR'YSTAR); C = chol(V); accept = 0; do while accept == 0; PHI_G = PHI + C'rndn(2,1); coef = -rev(PHI_G)|1; root = polyroot(coef); rootmod = abs(root); if minc(rootmod) ge 1.0001; accept = 1; else; counter = counter + 1; accept = 0; endif; endo; YSTAR=ZTT[5:t,1]- PHI_G[1,1]*ZTT[4:T-1] - PHI_G[2,1]*ZTT[3:T-2]; XSTAR=ONES(T-4,1)~(STT[5:T,1]-PHI_G[1]*STT[4:T-1]-PHI_G[2]*STT[3:T-2]); V = invpd(R0_m + XSTAR'XSTAR); MU = V*(R0_*t0_m + XSTAR'YSTAR); C = chol(V); accept = 0; do while accept == 0; MU_G = MU + C'rndn(2,1); IF MU_G[2] GT 0; accept = 1; endif; endo; RETP(PHI_G[1],PHI_G[2],(MU_G[1]/(1-phi_g[1]-phi_g[2])),MU_G[2]); ENDP; @=========================================================================@ @=======A PROCEDURE THAT GENERATES Zt ===================================@ @=========================================================================@ PROC GEN_ZT; local ppr,qpr,phi1,phi2,mu_0,mu_1,sig_v,alpha1,psi1_1,psi1_2,sig_1,alpha2, psi2_1,psi2_2,sig_2,alpha3,psi3_1,psi3_2,sig_3,alpha4,psi4_1,psi4_2, sig_4,aaa,h,qq,rr,mm0,mm1,xtt_m,ptt_m,m_gen_m,v_gen_m,prob__1,prob__0, p_cf_0,p_cf_1,p_cf,t_vr,p_vr,pst_cf,pst_vr,f_cast,ss,k_gn,j_iter,jj, f_cast0,ss0,z_mat,x_tmp,p_tmp,g,ztstar,zt,pr_val,k_gn0,YYSTAR,tstar, st,alpha41,alpha42,alpha43,rnd_mat,JJJ; @For common component@ PPR=1-B1TT; @Pr[St=1/St-1=1]@ QPR=1-A1TT; @Pr[St=0/St-1=0]@ phi1=phi1tt; phi2=phi2tt; mu_0=mu0tt; mu_1=mu1tt; sig_v=1; @For 1st variable@ alpha1=lamd1tt; psi1_1=psi11tt; psi1_2=psi12tt; sig_1=sig21tt; @For 2nd variable@ alpha2=lamd2tt; psi2_1=psi21tt; psi2_2=psi22tt; sig_2=sig22tt; @For 3rd variable@ alpha3=lamd3tt; psi3_1=psi31tt; psi3_2=psi32tt; sig_3=sig23tt; @For 4th variable@ alpha4=lamd4tt; psi4_1=psi41tt; psi4_2=psi42tt; sig_4=sig24tt; ALPHA41=LAM41TT;ALPHA42=LAM42TT; ALPHA43=LAM43TT; @---------------------------------------------------@ AAA=(phi1~phi2~0~0~0~0)|(EYE(5)~ZEROS(5,1)); H= (alpha1~-alpha1*psi1_1~-alpha1*psi1_2~0~0~0)| (alpha2~-alpha2*psi2_1~-alpha2*psi2_2~0~0~0)| (alpha3~-alpha3*psi3_1~-alpha3*psi3_2~0~0~0)| (alpha4~ -ALPHA4*PSI4_1+ALPHA41~ -ALPHA4*PSI4_2-ALPHA41*PSI4_1+ALPHA42~ -ALPHA41*PSI4_2-ALPHA42*PSI4_1+ALPHA43~ -ALPHA42*PSI4_2-ALPHA43*PSI4_1~ -ALPHA43*PSI4_2); QQ=zeros(6,6); QQ[1,1]=sig_v; RR=zeros(4,4); RR[1,1]=SIG_1; RR[2,2]=SIG_2; RR[3,3]=SIG_3; RR[4,4]=SIG_4; YYSTAR=ZEROS(4,T-2); YYSTAR[1,.] = YY[1,3:T] -PSI1_1*YY[1,2:T-1] - PSI1_2*YY[1,1:T-2]; YYSTAR[2,.] = YY[2,3:T] -PSI2_1*YY[2,2:T-1] - PSI2_2*YY[2,1:T-2]; YYSTAR[3,.] = YY[3,3:T] -PSI3_1*YY[3,2:T-1] - PSI3_2*YY[3,1:T-2]; YYSTAR[4,.] = YY[4,3:T] -PSI4_1*YY[4,2:T-1] - PSI4_2*YY[4,1:T-2]; ST=STT[3:T,1]; TSTAR=T-2; mm0=mu_0|ZEROS(5,1); mm1=mu_1|ZEROS(5,1); xtt_m=zeros(6,TSTAR); ptt_m=zeros(6,(TSTAR)*6); @>>>>>>>>>>>>>>>>>>>>>>>>> INITIAL PROB. Pr[S0/Y0] @ PROB__1=(1-QPR)/(2-PPR-QPR); @Pr[St-1=1/Yt-1], STEADY STATE PROB.@ PROB__0=1-PROB__1 ; @Pr[ST-1=0/Yt-1], STEADY STATE PROB @ p_cf=(mu_0+mu_1*prob__1)*ones(6,1); t_vr=inv(eye(6^2)-AAA.*.AAA)*VEC(QQ); p_vr=t_vr[1:6,1]~t_vr[7:12,1]~t_vr[13:18,1]~ t_vr[19:24,1]~t_vr[25:30,1]~t_vr[31:36,1]; J_ITER = 3; DO UNTIL J_ITER>TSTAR; @===================FILTERING=============================@ PST_CF = AAA * P_CF + (mm0 + mm1*St[j_iter,1]) -phi1*(mm0+mm1*st[j_iter-1,1]) -phi2*(mm0+mm1*st[j_iter-2,1]); PST_VR = AAA * P_VR * AAA' +QQ; F_CAST= YYSTAR[.,J_ITER]- H * PST_CF; SS = H * PST_VR * H' + RR; K_GN =PST_VR * H'*inv(SS); P_CF = PST_CF + K_GN * F_CAST; P_VR = (EYE(6) - K_GN * H ) * PST_VR; @=========================================================@ XTT_M[.,J_ITER]=P_CF; PTT_M[.,(J_ITER*6-5):(J_ITER*6)]=P_VR; J_ITER = J_ITER+1; ENDO; @-----------------------------------------------------------------@ Z_MAT=ZEROS(TSTAR,1); RND_MAT=RNDN(TSTAR+10,1); X_TMP = XTT_M[.,TSTAR]; P_TMP = PTT_M[.,(TSTAR*6-5):TSTAR*6]; G = CHOL(P_TMP)'; @P_TMP=G*G'@ ZTSTAR = RND_MAT[1:6,1]; ZT = X_TMP + G*ZTSTAR; Z_MAT[TSTAR,1] = ZT[1,1]; J_ITER=TSTAR-1; DO UNTIL J_ITER<3; JJJ=TSTAR-J_ITER+1; JJ=1; DO UNTIL JJ>1;@ROWS(P_CF);@ F_CAST0 = ZT[JJ,1] - AAA[JJ,.]*XTT_M[.,J_ITER] -( (mu_0 + mu_1*St[j_iter+1,1]) -phi1*(mu_0+mu_1*st[j_iter-1+1,1]) -phi2*(mu_0+mu_1*st[j_iter-2+1,1])); SS0 = AAA[JJ,.]*PTT_M[.,(J_ITER*6-5):J_ITER*6] *AAA[JJ,.]' + QQ[JJ,JJ]; K_GN0 =PTT_M[.,(J_ITER*6-5):J_ITER*6] * AAA[JJ,.]'/SS0; XTT_M[.,J_ITER] = XTT_M[.,J_ITER] + K_GN0*F_CAST0; PTT_M[.,(J_ITER*6-5):J_ITER*6] = (EYE(6) - K_GN0*AAA[JJ,.])*PTT_M[.,(J_ITER*6-5):J_ITER*6]; JJ=JJ+1; ENDO; X_TMP = XTT_M[.,J_ITER]; P_TMP = PTT_M[.,(J_ITER*6-5):J_ITER*6]; G = CHOL(P_TMP)'; @P_TMP=G*G'@ ZTSTAR =RND_MAT[JJJ:JJJ+5,1]; ZT = X_TMP + G*ZTSTAR; Z_MAT[J_ITER,1]=ZT[1,1]; J_ITER=J_ITER-1; ENDO; Z_MAT=REV(ZT[2:5,1])|Z_MAT[3:tstar,1]; RETP(Z_MAT); ENDP; @=========================================================================@ @====PROCEDURE THAT GENERATES ST==========================================@ @=========================================================================@ PROC GEN_ST; local lag_ar,no_st,dmnsion,st_mat,j,st1,st,t0,ystar,x_mat,tstar,flt_pr, ppr,qpr,phi1,phi2,mu_0,mu_1,sig_v,st2, var_l,pr_trf0,pr_trf,f1,mu_mat,prob__1,prob__0,prob__t,prob__, j_iter,f_cast1,prob_dd,pr_vl,pr_val,pro_,s_t,p0,p1; LAG_AR=2; NO_ST=lag_ar+1; @ NUMBER OF STATES TO BE CONSIDERED@ DMNSION=2^NO_ST; st_mat=zeros(DMNSION,NO_ST); j=1; st2=0;do until st2>1; st1=0; do until st1>1; @ S_{t-1} S_t @ st=0; do until st>1; @ 0 0 @ @ 0 1 @ st_mat[j,.]=st2~st1~st; @ 1 0 @ @ 1 1 @ j=j+1; st=st+1;endo; st1=st1+1;endo; st2=st2+1;endo; T0=ROWS(ZTT); YSTAR=ZTT[5:T0,1]; X_MAT=ZTT[4:T0-1,1]~ZTT[3:T0-2]; TSTAR=ROWS(YSTAR); FLT_PR=ZEROS(2,TSTAR); @For common component@ PPR=1-B1TT; @Pr[St=1/St-1=1]@ QPR=1-A1TT; @Pr[St=0/St-1=0]@ phi1=phi1tt; phi2=phi2tt; mu_0=mu0tt; mu_1=mu1tt; sig_v=1; VAR_L=SIG_V*ONES(DMNSION,1); pr_trf0=qpr|(1-qpr)|(1-ppr)|ppr; @transition prob. mat.@ pr_trf=pr_trf0|pr_trf0; @2^3x1@ @<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@ F1=YSTAR - X_MAT*(PHI1|PHI2); mu_mat=mu_0*ones(dmnsion,no_st) + st_mat*mu_1; @++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++@ @BEGIN FILTERING ++++++++++++++++++++++++++++++++++++++++++++++++@ @++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++@ PROB__1=(1-QPR)/(2-PPR-QPR); @Pr[St-1=1/Yt-1], STEADY STATE PROB.@ PROB__0=1-PROB__1 ; @Pr[ST-1=0/Yt-1], STEADY STATE PROB @ PROB__T=prob__0|prob__1; @4x1@ PROB__T=VECR(PROB__T~PROB__T).*PR_TRF0; prob__=VECR(PROB__T~PROB__T); J_ITER=1; DO UNTIL J_ITER>TSTAR; f_cast1=f1[j_iter,1]*ones(DMNSION,1)- mu_mat*(-phi2|-phi1|1); prob_dd=PR_TRF .* PROB__; pr_vl=(1./SQRT(2.*PI.*var_L)).*EXP(-0.5*f_cast1.*f_cast1./var_L) .*prob_dd; @PR[St-m,...,St-1,St|Y_t]@ pr_val=sumc(pr_vl); pro_=pr_vl/pr_val; PROB__T=PRO_[1:DMNSION/2,1]+PRO_[DMNSION/2+1:DMNSION,1]; @Pr[St-m+1,...,St/Yt]@ PROB__=VECR(PROB__T~PROB__T); PROB__T=PROB__T[1:ROWS(PROB__T)/2,1]+ PROB__T[ROWS(PROB__T)/2+1:ROWS(PROB__T),1]; FLT_PR[.,J_ITER]=PROB__T; @Pr[St|Yt], 2x1@ J_ITER = J_ITER+1; ENDO; @==========GENERATE S_T BASED ON CARTER & KOHN (1994)===========@ S_T=ZEROS(TSTAR,1); S_T[TSTAR,1]= BINGEN(FLT_PR[1,TSTAR],FLT_PR[2,TSTAR],1); J_ITER=TSTAR-1; DO UNTIL J_ITER<1; IF S_T[J_ITER+1,1]==0; P0=QPR*FLT_PR[1,J_ITER]; P1=(1-PPR)*FLT_PR[2,J_ITER]; ELSEIF S_T[J_ITER+1,1]==1; P0=(1-QPR)*FLT_PR[1,J_ITER]; P1=PPR*FLT_PR[2,J_ITER]; ENDIF; S_T[J_ITER,1]=BINGEN(P0,P1,1); J_ITER=J_ITER-1; ENDO; RETP(1|1|1|1|S_T); ENDP; @++++++++++++++++++++++++++++++++++++++++++++++++++++++++@ @=========================================================================@ @=======A PROCEDURE THAT CALCULATES STEADY-STATE KALMAN GAIN,=============@ @=======GIVEN S_t {t=1,,T}, AND OTHER PARAMETERS OF THE MODEL=============@ PROC SS_KGG(tmpry); local ppr,qpr,phi1,phi2,mu_0,mu_1,sig_v,alpha1,psi1_1,psi1_2,sig_1,alpha2, psi2_1,psi2_2,sig_2,alpha3,psi3_1,psi3_2,sig_3,alpha4,psi4_1,psi4_2, sig_4,aaa,h,qq,rr,mm0,mm1,xtt_m,ptt_m,m_gen_m,v_gen_m,prob__1,prob__0, p_cf_0,p_cf_1,p_cf,t_vr,p_vr,pst_cf,pst_vr,f_cast,ss,k_gn,j_iter,jj, f_cast0,ss0,z_mat,x_tmp,p_tmp,g,ztstar,zt,pr_val,k_gn0,YSTAR,tstar, st,tmpm0,tmpm1,alpha41,alpha42,alpha43; @For common component@ PPR=1-B1TT; @Pr[St=1/St-1=1]@ QPR=1-A1TT; @Pr[St=0/St-1=0]@ phi1=phi1tt; phi2=phi2tt; mu_0=mu0tt; mu_1=mu1tt; sig_v=1; @For 1st variable@ alpha1=lamd1tt; psi1_1=psi11tt; psi1_2=psi12tt; sig_1=sig21tt; @For 2nd variable@ alpha2=lamd2tt; psi2_1=psi21tt; psi2_2=psi22tt; sig_2=sig22tt; @For 3rd variable@ alpha3=lamd3tt; psi3_1=psi31tt; psi3_2=psi32tt; sig_3=sig23tt; @For 4th variable@ alpha4=lamd4tt; psi4_1=psi41tt; psi4_2=psi42tt; sig_4=sig24tt; alpha41=lam41tt;alpha42=lam42tt; alpha43=lam43tt; @---------------------------------------------------@ AAA=(phi1~phi2~0~0~zeros(1,8))| (1~ 0~0~0~zeros(1,8))| (0~ 1~0~0~zeros(1,8))| (0~ 0~1~0~zeros(1,8))| (zeros(1,4)~psi1_1~psi1_2~zeros(1,6))| (zeros(1,4)~1~zeros(1,7))| (zeros(1,6)~psi2_1~psi2_2~zeros(1,4))| (zeros(1,6)~1~zeros(1,5))| (zeros(1,8)~psi3_1~psi3_2~zeros(1,2))| (zeros(1,8)~1~zeros(1,3))| (zeros(1,10)~psi4_1~psi4_2)| (zeros(1,10)~1~zeros(1,1)); H= (alpha1~zeros(1,3)~1~zeros(1,7))| (alpha2~zeros(1,5)~1~zeros(1,5))| (alpha3~zeros(1,7)~1~zeros(1,3))| (alpha4~alpha41~alpha42~alpha43~zeros(1,6)~1~zeros(1,1)); QQ= zeros(12,12); QQ[1,1]=sig_v; QQ[5,5]=sig_1; QQ[7,7]=sig_2; QQ[9,9]=sig_3; QQ[11,11]=sig_4; RR=ZEROS(4,4); ST=STT[1:T,1]; TSTAR=T; YSTAR=YY; PROB__1=(1-QPR)/(2-PPR-QPR); @Pr[St-1=1/Yt-1], STEADY STATE PROB.@ PROB__0=1-PROB__1 ; @Pr[ST-1=0/Yt-1], STEADY STATE PROB @ mm0=mu_0|zeros(11,1); mm1=mu_1|zeros(11,1); p_cf=(mu_0+mu_1*prob__1)*ones(4,1)|zeros(8,1); t_vr=inv(eye(12^2)-AAA.*.AAA)*VEC(QQ); p_vr=t_vr[1:12,1]~t_vr[13:24,1]~t_vr[25:36,1]~ t_vr[37:48,1]~t_vr[49:60,1]~t_vr[61:72,1]~ t_vr[73:84,1]~t_vr[85:96,1]~t_vr[97:108,1]~ t_vr[109:120,1]~t_vr[121:132,1]~t_vr[133:144,1]; J_ITER = 3; DO UNTIL J_ITER>TSTAR; @===================FILTERING=============================@ PST_CF = AAA * P_CF + (mm0 + mm1*St[j_iter,1]) -phi1*(mm0+mm1*st[j_iter-1,1])-phi2*(mm0+mm1*st[j_iter-2,1]); PST_VR = AAA * P_VR * AAA' +QQ; F_CAST= YSTAR[.,J_ITER]- H * PST_CF; SS = H * PST_VR * H' + RR; K_GN =PST_VR * H'*inv(SS); P_CF = PST_CF + K_GN * F_CAST; P_VR = (EYE(12) - K_GN * H ) * PST_VR; @=========================================================@ J_ITER = J_ITER+1; ENDO; RETP(K_GN); ENDP; @=========================================================================@ @=======A PROCEDURE THAT ESTIMATES COMPOSITE COINCIDENT INDEX, ===========@ @=======GIVEN S_t {t=1,2,...,T} AND OTHER SIMULATED PARAMETERS============@ PROC (2)= GEN_CCI; local ppr,qpr,phi1,phi2,mu_0,mu_1,sig_v,alpha1,psi1_1,psi1_2,sig_1,alpha2, psi2_1,psi2_2,sig_2,alpha3,psi3_1,psi3_2,sig_3,alpha4,psi4_1,psi4_2, sig_4,aaa,h,qq,rr,mm0,mm1,xtt_m,ptt_m,m_gen_m,v_gen_m,prob__1,prob__0, p_cf_0,p_cf_1,p_cf,t_vr,p_vr,pst_cf,pst_vr,f_cast,ss,k_gn,j_iter,jj, f_cast0,ss0,z_mat,x_tmp,p_tmp,g,ztstar,zt,pr_val,k_gn0,YYSTAR,tstar, st,a_tmp,q_tmp,h_tmp,t_vr0,p_vr0,ss_kg,kk,alpha,c_bar,delta,delta_m, beta,CTT_MAT,yy_m0,yy_m1,alp0,alp1,c_bar0,c_barM,beta0,beta1,ALPHA41, ALPHA42,ALPHA43; @For common component@ PPR=1-B1TT; @Pr[St=1/St-1=1]@ QPR=1-A1TT; @Pr[St=0/St-1=0]@ phi1=phi1tt; phi2=phi2tt; mu_0=mu0tt; mu_1=mu1tt; sig_v=1; @For 1st variable@ alpha1=lamd1tt; psi1_1=psi11tt; psi1_2=psi12tt; sig_1=sig21tt; @For 2nd variable@ alpha2=lamd2tt; psi2_1=psi21tt; psi2_2=psi22tt; sig_2=sig22tt; @For 3rd variable@ alpha3=lamd3tt; psi3_1=psi31tt; psi3_2=psi32tt; sig_3=sig23tt; @For 4th variable@ alpha4=lamd4tt; psi4_1=psi41tt; psi4_2=psi42tt; sig_4=sig24tt; alpha41=lam41tt;alpha42=lam42tt; alpha43=lam43tt; @---------------------------------------------------@ AAA=(phi1~phi2~0~0~zeros(1,8))| (1~ 0~0~0~zeros(1,8))| (0~ 1~0~0~zeros(1,8))| (0~ 0~1~0~zeros(1,8))| (zeros(1,4)~psi1_1~psi1_2~zeros(1,6))| (zeros(1,4)~1~zeros(1,7))| (zeros(1,6)~psi2_1~psi2_2~zeros(1,4))| (zeros(1,6)~1~zeros(1,5))| (zeros(1,8)~psi3_1~psi3_2~zeros(1,2))| (zeros(1,8)~1~zeros(1,3))| (zeros(1,10)~psi4_1~psi4_2)| (zeros(1,10)~1~zeros(1,1)); H= (alpha1~zeros(1,3)~1~zeros(1,7))| (alpha2~zeros(1,5)~1~zeros(1,5))| (alpha3~zeros(1,7)~1~zeros(1,3))| (alpha4~alpha41~alpha42~alpha43~zeros(1,6)~1~zeros(1,1)); QQ= zeros(12,12); QQ[1,1]=sig_v; QQ[5,5]=sig_1; QQ[7,7]=sig_2; QQ[9,9]=sig_3; QQ[11,11]=sig_4; RR=ZEROS(4,4); ST=STT[1:T,1]; TSTAR=T; @>>>>>>>>>>>>>>>>>>>>>>>>> INITIAL PROB. Pr[S0/Y0] @ PROB__1=(1-QPR)/(2-PPR-QPR); @Pr[St-1=1/Yt-1], STEADY STATE PROB.@ PROB__0=1-PROB__1 ; @Pr[ST-1=0/Yt-1], STEADY STATE PROB @ {SS_KG}=SS_KGG(100); KK=(EYE(12)-SS_KG*H)*AAA; ALP0=INV(EYE(12)-KK)*SS_KG*(YY_M); C_BAR= ALP0[1,1]; C_BARM=C_BAR|ZEROS(10,1); c_bar; ctt_mat=zeros(t,1); j_iter=3; do until j_iter>tstar; ctt_mat[j_iter,1]=ztt[j_iter,1]+ctt_mat[j_iter-1,1]+c_bar; @C_t|T@ j_iter=j_iter+1; endo; /* BETA= YY_M - (ALPHA1|ALPHA2|ALPHA3|ALPHA4)*C_BAR; p_cf=(mu_0+mu_1*prob__1+c_bar)*ones(2,1)|ZEROS(8,1)|0; mm0=(mu_0)|ZEROS(10,1); mm1=(mu_1)|ZEROS(10,1); CTT_MAT=ZEROS(T,1); J_ITER = 3; DO UNTIL J_ITER>TSTAR; @===================FILTERING=============================@ PST_CF = AAA * P_CF + (mm0 + mm1*St[j_iter,1]+c_barm) -phi1*(mm0+mm1*st[j_iter-1,1]+c_barm) -phi2*(mm0+mm1*st[j_iter-2,1]+c_barm); PST_VR = AAA * P_VR * AAA' +QQ; F_CAST= YYY[.,J_ITER]- H * PST_CF -beta; SS = H * PST_VR * H' + RR; K_GN =PST_VR * H'*inv(SS); P_CF = PST_CF + K_GN * F_CAST; P_VR = (EYE(11) - K_GN * H ) * PST_VR; @=========================================================@ CTT_MAT[J_ITER,1]=P_CF[1,1]+P_CF[11,1]; @C_t|t@ J_ITER = J_ITER+1; ENDO; */ RETP(CTT_MAT, c_bar); ENDP; @=============================================================@ proc sw_cnt(s); local n,m,switch,tt,st1,st,swch; swch=zeros(t,4); tt = 2; do until tt>t; if s[tt-1,1]==0 and s[tt,1]==0; swch[tt,1]=1; elseif s[tt-1,1]==0 and s[tt,1]==1; swch[tt,2]=1; elseif s[tt-1,1]==1 and s[tt,1]==0; swch[tt,3]=1; elseif s[tt-1,1]==1 and s[tt,1]==1; swch[tt,4]=1; endif; tt = tt + 1; endo; retp(swch); endp; @=============================================================@ @========================================================================@ @========================================================================@ /* ADDPROCS.GSI */ /* compiled file used in Albert and Chib (JBES, 1993, 11,1-15) */ /* please contact Siddhartha Chib at chib@olin.wustl.edu if there are problems*/ proc rndc(a); /* chisquare */ local x,w; a = a/2; w= rg1(a); x = w * 2; retp(x); endp; proc rndb(a,b); /* beta */ local x,a1n,a1d; a1n = rg1(a); a1d = rg1(b); x = a1n / (a1n + a1d); retp(x); endp; proc rg1(a); local x,j,w,u; if a gt 1; x = rg2(a); elseif a lt 1; a = a + 1; u = rndu(1,1); x = rg2(a)*u^(1/a); elseif a == 1; x = -ln(rndu(1,1)); endif; retp(x); endp; proc rg2(a); local gam,accept,b,c,j,x,z,u,v,w,y; b = a-1; c = 3*a - .75; accept = 0; do while accept == 0; u = rndu(1,1); v = rndu(1,1); w = u*(1-u); y = sqrt(c/w)*(u-.5); x = b+y; if x ge 0; z = 64*(w^3)*(v^2); accept = z le ( 1-(2*y^2)/x ); if accept == 0; accept = ln(z) le 2*(b*ln(x/b) - y); endif; endif; endo; retp(x); endp; @================================================@ /* proc gendiff(z,rho); local ztrim,i,zgdiff; i = 1; ztrim = trimr(z,r,0); zgdiff = ztrim; do while i le r; zgdiff = zgdiff - rho[i]*trimr(z,r-i,i); i = i + 1; endo; retp(zgdiff); endp; */ proc pdfn1(yt,m,s2); local z,p; z = (yt - m)/sqrt(s2); p = pdfn(z)/sqrt(s2); retp(p); endp; trace 0; proc bingen(p0,p1,m); local pr0,s,u; pr0 = p0/(p0+p1); /* prob(s=0) */ u = rndu(m,1); s = u .ge pr0; retp(s); endp; proc markov(a,b,ss0,n); local s,u,i ; s = zeros(n,1); s[1,1] = ss0; i = 2; do while i le n; if s[i-1,1] == 0; u = rndu(1,1); if u le a ; s[i,1] = 1; else; s[i,1] = 0; endif; elseif s[i-1,1] == 1; u = rndu(1,1); if u le (1-b) ; s[i,1] = 1; else; s[i,1] = 0; endif; endif; i = i + 1; endo; retp(s); endp; proc switchg(s,g); local n,m,switch,t,st1,st; n = rows(s); m = rows(g); switch = zeros(m,m); /* to store the transitions */ t = 2; do while t le n; st1 = s[t-1]; st = s[t]; switch[st1,st] = switch[st1,st] + 1; t = t+ 1; endo; retp(switch); endp;