/* =============================================================== File: a:\srfms.opt MODEL: Fixed Parameter model with Markov Switching Variance and Mean for stock return predictability. r(t,t+1) - MU_t = beta_k * (r(t-k,t) - k * MU_t-j) + e_t e_t ~ N(0, RR_t) RR_t = SIG0^2(1-S_t) + SIG1^2*S_t MU_t = MU0 + MU1*S_t Pr[S_t=1|S_t-1=1]=p, Pr[S_t=0|S_t-1=0]=q ===> 6 parameters (MU0, MU1, beta_k, sig0, sig1) to be estimated WRITTEN BY JAMES MORLEY, DEPT. OF ECONOMICS, UNIVERSITY OF WASHINGTON MODIFIED FROM A PROGRAM BY CHANG_JIN KIM, DEPT. OF ECONOMICS, KOREA UNIVERSITY ================================================================ */ @ Version 5/8/98@ New; library optmum,pgraph; #include optmum.ext; optset; format /m1 /rd 8,6; k=120; @Select holding period (in months) for k-period returns@ r=2; @Select equal (r=1) or value (r=2) weighted returns@ @Select decile(j-2) (r=j)@ load data[852,12]=c:\CRSPD96.DAT; /* column 1: CRSP monthly index of NYSE equal-weighted total real return 2: CRSP monthly index of NYSE value-weighted total real return 3-12: decile returns Sample period -- 1926:1--1986:12,... 732 observations */ T=rows(data); @Create the Data@ @----------------------------------------------------------------@ @Equal Weighted one month and lagged k-period returns@ ew=ln(1+data[1:T,1]); ew_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; ew_kl[j+1,1]=sumc(ew[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Value Weighted one month and lagged k-period returns@ vw=ln(1+data[1:T,2]); vw_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; vw_kl[j+1,1]=sumc(vw[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 1 one month and lagged k-period returns@ d1=ln(1+data[1:T,3]); d1_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d1_kl[j+1,1]=sumc(d1[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 2 one month and lagged k-period returns@ d2=ln(1+data[1:T,4]); d2_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d2_kl[j+1,1]=sumc(d2[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 3 one month and lagged k-period returns@ d3=ln(1+data[1:T,5]); d3_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d3_kl[j+1,1]=sumc(d3[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 4 one month and lagged k-period returns@ d4=ln(1+data[1:T,6]); d4_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d4_kl[j+1,1]=sumc(d4[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 5 one month and lagged k-period returns@ d5=ln(1+data[1:T,7]); d5_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d5_kl[j+1,1]=sumc(d5[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 6 one month and lagged k-period returns@ d6=ln(1+data[1:T,8]); d6_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d6_kl[j+1,1]=sumc(d6[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 7 one month and lagged k-period returns@ d7=ln(1+data[1:T,9]); d7_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d7_kl[j+1,1]=sumc(d7[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 8 one month and lagged k-period returns@ d8=ln(1+data[1:T,10]); d8_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d8_kl[j+1,1]=sumc(d8[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 9 one month and lagged k-period returns@ d9=ln(1+data[1:T,11]); d9_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d9_kl[j+1,1]=sumc(d9[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 10 one month and lagged k-period returns@ d10=ln(1+data[1:T,12]); d10_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d10_kl[j+1,1]=sumc(d10[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @------------------------------------------------------------------@ output file=c:\p6\r2k120ms.out reset; @S=1->1926-86,S=2->1926-46,S=3->1947-86@ output off; data=ew~vw~d1~d2~d3~d4~d5~d6~d7~d8~d9~d10~ew_kl~vw_kl~ d1_kl~d2_kl~d3_kl~d4_kl~d5_kl~d6_kl~d7_kl~d8_kl~d9_kl~d10_kl; data=data[k+1:T,1:24]; S=T-k; @Subsample after adjusting for lagged X@ if r<2; yy=data[.,1]; xx=data[.,13]; @ew returns@ goto skip1; endif; if r<3; yy=data[.,2]; xx=data[.,14]; @vw returns@ goto skip1; endif; if r<4; yy=data[.,3]; xx=data[.,15]; @vw returns@ goto skip1; endif; if r<5; yy=data[.,4]; xx=data[.,16]; @vw returns@ goto skip1; endif; if r<6; yy=data[.,5]; xx=data[.,17]; @vw returns@ goto skip1; endif; if r<7; yy=data[.,6]; xx=data[.,18]; @vw returns@ goto skip1; endif; if r<8; yy=data[.,7]; xx=data[.,19]; @vw returns@ goto skip1; endif; if r<9; yy=data[.,8]; xx=data[.,20]; @vw returns@ goto skip1; endif; if r<10; yy=data[.,9]; xx=data[.,21]; @vw returns@ goto skip1; endif; if r<11; yy=data[.,10]; xx=data[.,22]; @vw returns@ goto skip1; endif; if r<12; yy=data[.,11]; xx=data[.,23]; @vw returns@ goto skip1; endif; if r<13; yy=data[.,12]; xx=data[.,24]; @vw returns@ goto skip1; endif; skip1: @================= Initialize Global Variables============@ START=1; PRMTR_NO=4*(k-1)+8; @ >>Number of parameters to be estimated. >>@ pr_intl={3.395345 2.588071 5.582229 -3.382499 3.905454 0.000713 -3.809880 -2.068871 }; pr_1=3.359121 ;pr_2=2.582513 ;pr_3=5.711422 ; pr_4=-3.006611 ;pr_5=3.346763 ;pr_6=0.000991 ; pr_7=-3.958776 ;pr_8=-2.233841 ; PRMTR_SI= {0.071282 -0.002135 0.937767 -0.003440 -0.032673 -0.002146 0.944895 -0.003439 -0.032543 -0.001673 0.962353 -0.003417 -0.032529 -0.001689 0.969207 -0.003431 -0.032566 -0.001836 1.053530 -0.003324 -0.032618 -0.001841 1.058372 -0.003323 -0.032668 -0.001813 1.063010 -0.003321 -0.032627 -0.001801 1.067416 -0.003294 -0.032703 -0.001819 1.071595 -0.003290 -0.032673 -0.001808 1.075536 -0.003288 -0.032688 -0.001823 1.079227 -0.003291 -0.032692 -0.001850 1.082691 -0.003270 -0.032751 -0.001854 1.085874 -0.003264 -0.032755 -0.001876 1.088785 -0.003264 -0.032761 -0.001875 1.091440 -0.003258 -0.032764 -0.001850 1.093807 -0.003238 -0.032996 -0.001765 1.069639 -0.002949 -0.033035 -0.001750 1.067188 -0.002947 -0.033093 -0.001755 1.064545 -0.002935 -0.033057 -0.001755 1.061715 -0.002921 -0.033108 -0.001726 1.058679 -0.002890 -0.033124 -0.001759 1.055475 -0.002886 -0.033141 -0.001739 1.052074 -0.002868 -0.033148 -0.001732 1.048533 -0.002862 -0.033129 -0.001706 1.044809 -0.002847 -0.033148 -0.001720 1.040933 -0.002824 -0.033141 -0.001726 1.036916 -0.002811 -0.033182 -0.001745 1.032764 -0.002779 -0.033169 -0.001710 1.028481 -0.002766 -0.033173 -0.001708 1.024062 -0.002748 -0.033198 -0.001691 1.019522 -0.002743 -0.033216 -0.001675 1.014884 -0.002735 -0.033223 -0.001674 1.010150 -0.002701 -0.033226 -0.001692 1.005317 -0.002684 -0.033205 -0.001683 1.000390 -0.002680 -0.033179 -0.001686 0.995394 -0.002674 -0.033194 -0.001682 0.990306 -0.002643 -0.033192 -0.001700 0.985175 -0.002639 -0.032761 -0.001875 1.091440 -0.003258 -0.032764 -0.001850 1.093807 -0.003238 -0.032996 -0.001765 1.069639 -0.002949 -0.033035 -0.001750 1.067188 -0.002947 -0.033093 -0.001755 1.064545 -0.002935 -0.033057 -0.001755 1.061715 -0.002921 -0.033108 -0.001726 1.058679 -0.002890 -0.033124 -0.001759 1.055475 -0.002886 -0.033141 -0.001739 1.052074 -0.002868 -0.033148 -0.001732 1.048533 -0.002862 -0.033129 -0.001706 1.044809 -0.002847 -0.033148 -0.001720 1.040933 -0.002824 -0.033141 -0.001726 1.036916 -0.002811 -0.033182 -0.001745 1.032764 -0.002779 -0.033169 -0.001710 1.028481 -0.002766 -0.033173 -0.001708 1.024062 -0.002748 -0.033198 -0.001691 1.019522 -0.002743 -0.033216 -0.001675 1.014884 -0.002735 -0.033223 -0.001674 1.010150 -0.002701 -0.033226 -0.001692 1.005317 -0.002684 -0.033205 -0.001683 1.000390 -0.002680 -0.033179 -0.001686 0.995394 -0.002674 -0.033194 -0.001682 0.990306 -0.002643 -0.033192 -0.001700 0.985175 -0.002639 -0.032761 -0.001875 1.091440 -0.003258 -0.032764 -0.001850 1.093807 -0.003238 -0.032996 -0.001765 1.069639 -0.002949 -0.033035 -0.001750 1.067188 -0.002947 -0.033093 -0.001755 1.064545 -0.002935 -0.033057 -0.001755 1.061715 -0.002921 -0.033108 -0.001726 1.058679 -0.002890 -0.033124 -0.001759 1.055475 -0.002886 -0.033141 -0.001739 1.052074 -0.002868 -0.033148 -0.001732 1.048533 -0.002862 -0.033129 -0.001706 1.044809 -0.002847 -0.033148 -0.001720 1.040933 -0.002824 -0.033141 -0.001726 1.036916 -0.002811 -0.033182 -0.001745 1.032764 -0.002779 -0.033169 -0.001710 1.028481 -0.002766 -0.033173 -0.001708 1.024062 -0.002748 -0.033198 -0.001691 1.019522 -0.002743 -0.033216 -0.001675 1.014884 -0.002735 -0.033223 -0.001674 1.010150 -0.002701 -0.033226 -0.001692 1.005317 -0.002684 -0.033205 -0.001683 1.000390 -0.002680 -0.033179 -0.001686 0.995394 -0.002674 -0.033194 -0.001682 0.990306 -0.002643 -0.033192 -0.001700 0.985175 -0.002639 -0.032761 -0.001875 1.091440 -0.003258 -0.032764 -0.001850 1.093807 -0.003238 -0.032996 -0.001765 1.069639 -0.002949 -0.033035 -0.001750 1.067188 -0.002947 -0.033093 -0.001755 1.064545 -0.002935 -0.033057 -0.001755 1.061715 -0.002921 -0.033108 -0.001726 1.058679 -0.002890 -0.033124 -0.001759 1.055475 -0.002886 -0.033141 -0.001739 1.052074 -0.002868 -0.033148 -0.001732 1.048533 -0.002862 -0.033129 -0.001706 1.044809 -0.002847 -0.033148 -0.001720 1.040933 -0.002824 -0.033141 -0.001726 1.036916 -0.002811 -0.033182 -0.001745 1.032764 -0.002779 -0.033169 -0.001710 1.028481 -0.002766 -0.033173 -0.001708 1.024062 -0.002748 -0.033198 -0.001691 1.019522 -0.002743 -0.033216 -0.001675 1.014884 -0.002735 -0.033223 -0.001674 1.010150 -0.002701 -0.033226 -0.001692 1.005317 -0.002684 -0.033205 -0.001683 1.000390 -0.002680 -0.033179 -0.001686 0.995394 -0.002674 -0.033194 -0.001682 0.990306 -0.002643 -0.033192 -0.001700 0.985175 -0.002639 -0.033229 -0.001701 0.979970 -0.002630 -0.033202 -0.001687 0.974727 -0.002589 -0.033206 -0.001679 0.969435 -0.002575 -0.033205 -0.001680 0.964105 -0.002560 -0.033196 -0.001699 0.958741 -0.002533 -0.033174 -0.001696 0.953362 -0.002511 -0.033203 -0.001677 0.947975 -0.002497 -0.033172 -0.001693 0.942553 -0.002480 -0.033170 -0.027193 -0.047673 }; @PRMTR_IN=pr_1|pr_2|pr_3|pr_4|pr_5|pr_6|pr_7|pr_8;@ @PRMTR_IN=pr_intl'|0.2*ones(4*(k-1)-1,1);@ pr_int=pr_intl'; PRMTR_IN=pr_intl; PRMTR_IN[4]=pr_intl[4]+pr_intl[5]; PRMTR_IN[5]=pr_intl[4]; @ PRMTR_IN=pr_intl'|prmtr_si';@ @OUTPUT@ output on; if r<2; "Equal Weighted Real Returns"; goto skip2; endif; if r<3; "Value Weighted Real Returns"; goto skip2; endif; if r<4; "Decile(1) Real Returns"; goto skip2; endif; if r<5; "Decile(2) Real Returns"; goto skip2; endif; if r<6; "Decile(3) Real Returns"; goto skip2; endif; if r<7; "Decile(4) Real Returns"; goto skip2; endif; if r<8; "Decile(5) Real Returns"; goto skip2; endif; if r<9; "Decile(6) Real Returns"; goto skip2; endif; if r<10; "Decile(7) Real Returns"; goto skip2; endif; if r<11; "Decile(8) Real Returns"; goto skip2; endif; if r<12; "Decile(9) Real Returns"; goto skip2; endif; if r<13; "Decile(10) Real Returns"; goto skip2; endif; skip2: "(2-state Markov-switching in variance and mean)"; "==MODEL======================================================="; " r(t,t+1) - MU_t = beta_k * (r(t-k,t) - k * MU_t-j) + e_t"; " e_t ~ N(0, RR_t)"; ""; " RR_t = (SIG0 + SIG1*S_t)^2"; " MU_t = MU0 + MU1*S_t"; " Pr[S_t=1|S_t-1=1]=p, Pr[S_t=0|S_t-1=0]=q"; "==============================================================="; " "; "==Estimation Details==========================================="; "k=";;k; "Start=";;START;;" => t=1 @ January, ";;1926+(k+Start-1)/12; " (first 12 observations are used to start up Kalman filter)"; "End=";;1926+T/12; "==============================================================="; " "; output off; @ Maximum Likelihood Estimation @ @_opgtol=1e-4;@ {xout,fout,gout,cout}=optmum(&lik_fcn,PRMTR_IN); @ prmtr estimates, -log lik value, Grandient, Hessian@ "Calculating Hessian..... Please be patient!!!!"; hout0=hessp(&lik_fcn,xout); hout=inv(hout0); PRM_FNL=TRANSA(xout); @ Estimated coefficients, constrained@ grdn_fnl=gradfd(&TRANSA,xout); Hsn_fnl=grdn_fnl*hout*grdn_fnl'; SD_fnl=sqrt(diag(Hsn_fnl)); @Standard errors of the estimated coefficients@ output on; "==FINAL OUTPUT================================================="; "likelihood value is ";; -fout; "Estimated parameters are:"; " sig0 sig1 mu0 mu1 ddd beta q p"; prm_fnl[1:8]'; "Standard errors of parameters are:"; sd_fnl[1:8]'; "==============================================================="; @"(";;xout';;")";@ "(";;xout';;PRMTR_SI;;")"; " "; output off; {pr_tl,pr_tt}=filter(xout); pr_tall=smooth(pr_tl,pr_tt); output file=c:\p6\r2k120ms.dta reset; pr_tt~pr_tall; output off; begwind; window(1,1,0); xlabel("Pr{S_t=1|Psi_T}"); xindex=seqa(1,1,rows(pr_tall)); title("Smoothed Probability of High Variance State"$+ ":\l[2-state Markov-switching Mean and Variance]"); XY(xindex,pr_tall); endwind; END; @ END OF MAIN PROGRAM @ @========================================================================@ @========================================================================@ PROC LIK_FCN(PRMTR1); LOCAL PRMTR, q,p,SIG0,SIG1,MU0,MU1,STP,PROB_LL,PROB_TL,s_index, PROB_TT1,PROB_TT2,PROB_TT3,f_ys,likv,j_iter,yt,xt,f_y, beta,ggg,PROBS_TL,ddd,MU1STAR; EXTERNAL PROC TRANS, TRANSA, V_PROB; PRMTR=TRANS(PRMTR1|PRMTR_SI'); @PRMTR=TRANSA(PRMTR1);@ @PRMTR=TRANS(PRMTR1);@ SIG0=PRMTR[1]; SIG1=PRMTR[2]; MU0=PRMTR[3]*ones(k-1,1); MU1STAR=PRMTR[4]*ones(k-1,1); DDD=PRMTR[5]; MU1=DDD; GGG=MU1STAR-DDD; beta=PRMTR[6]; q=PRMTR[7]; p=PRMTR[8]; STP=(1-SUMC(PRMTR[8+1:8+(k-1)*4-1]))|PRMTR[8+1:8+(k-1)*4-1]; @STP=ones(4*(k-1),1)/(4*(k-1));@ LOCATE 19,1; PRMTR[1:8]'; /* HAMILTON FILTER */ /*==============================================================*/ s_index=cumsumc(ones(k-1,1))-ones(k-1,1); @S_t*=S_t-2 + S_t-3 + ... + S_t-k+1@ @Startup Probability for Hamilton filter@ @Pr[S_t-1,S_t-k,S_t*|Y_t-1]@ @ 00 01 10 11 by k-1 @ PROB_LL=reshape(stp,k-1,4); PROBS_TL=(1-q)/(2-p-q)*ones(T,1); likv=0.0; j_iter=1; do until j_iter>S; yt=yy[j_iter,.]*ones(k-1,1); xt=xx[j_iter,.]*ones(k-1,1); @CALCULATE PROBABILITIES GIVEN INFORMATION UP TO T-1@ @Pr[S_t,S_t-1,S_t-k,S_t*|Y_t-1]=Pr[S_t|S_t-1]*Pr[S_t-1,S_t-k,S_t*|Y_t-1]@ @000 001 010 011 100 101 110 111 by k-1@ PROB_TL= q*PROB_LL[.,1]~q*PROB_LL[.,2]~(1-p)*PROB_LL[.,3]~(1-p)*PROB_LL[.,4]~ (1-q)*PROB_LL[.,1]~(1-q)*PROB_LL[.,2]~p*PROB_LL[.,3]~p*PROB_LL[.,4]; @Pr[S_t=1|Y_t-1]@ PROBS_TL[k+j_iter,1]=sumc(vecr(PROB_TL[.,5:8])); @CALCULATE MARGINAL DENSITY CONDITIONAL ON STATES@ @f(y_t|S_t,S_t-1,S_t-k,S_t*,Y_t-1)@ f_ys=V_PROB(SIG0^2,yt-MU0-PROBS_TL[k+j_iter,1]*GGG- beta*(xt-k*MU0-s_index[.,1].*MU1-sumc(PROBS_TL[j_iter:k-1+j_iter,1])*GGG))~ V_PROB(SIG0^2,yt-MU0-PROBS_TL[k+j_iter,1]*GGG-beta* (xt-k*MU0-(s_index[.,1]+ones(k-1,1)).*MU1- sumc(PROBS_TL[j_iter:k-1+j_iter,1])*GGG))~ V_PROB(SIG0^2,yt-MU0-PROBS_TL[k+j_iter,1]*GGG-beta* (xt-k*MU0-(s_index[.,1]+ones(k-1,1)).*MU1- sumc(PROBS_TL[j_iter:k-1+j_iter,1])*GGG))~ V_PROB(SIG0^2,yt-MU0-PROBS_TL[k+j_iter,1]*GGG-beta* (xt-k*MU0-(s_index[.,1]+2*ones(k-1,1)).*MU1- sumc(PROBS_TL[j_iter:k-1+j_iter,1])*GGG))~ V_PROB(SIG1^2,yt-MU0-MU1-PROBS_TL[k+j_iter,1]*GGG-beta* (xt-k*MU0-s_index[.,1].*MU1-sumc(PROBS_TL[j_iter:k-1+j_iter,1])*GGG))~ V_PROB(SIG1^2,yt-MU0-MU1-PROBS_TL[k+j_iter,1]*GGG-beta* (xt-k*MU0-(s_index[.,1]+ones(k-1,1)).*MU1- sumc(PROBS_TL[j_iter:k-1+j_iter,1])*GGG))~ V_PROB(SIG1^2,yt-MU0-MU1-PROBS_TL[k+j_iter,1]*GGG-beta* (xt-k*MU0-(s_index[.,1]+ones(k-1,1)).*MU1- sumc(PROBS_TL[j_iter:k-1+j_iter,1])*GGG))~ V_PROB(SIG1^2,yt-MU0-MU1-PROBS_TL[k+j_iter,1]*GGG-beta* (xt-k*MU0-(s_index[.,1]+2*ones(k-1,1)).*MU1- sumc(PROBS_TL[j_iter:k-1+j_iter,1])*GGG)); @CALCULATE MARGINAL DENSITY BY SUMMING ACROSS STATES@ @f(y_t|Y_t-1)@ f_y=sumc(vecr(f_ys.*PROB_TL)); @CALCULATE PROBABILITIES GIVEN INFORMATION UP TO T@ @Pr[S_t,S_t-1,S_t-k,S_t*|Y_t]@ @000 001 010 011 100 101 110 111 by k-1@ PROB_TT1=(f_ys.*PROB_TL)/f_y; @CALCULATE MORE PROBABILITIES GIVEN INFORMATION UP TO T@ @Pr[S_t,S_t-1,S_t-k,S_t-k+1,S_t*|Y_t]= Pr[S_t-k+1|S_t-k]*Pr[S_t,S_t-1,S_t-k,S_t*|Y_t]@ @0000 0010 0100 0110 1000 1010 1100 1110 0001 0011 0101 0111 1001 1011 1101 1111 by k-1@ PROB_TT2= q*PROB_TT1[.,1]~(1-p)*PROB_TT1[.,2]~q*PROB_TT1[.,3]~(1-p)*PROB_TT1[.,4]~ q*PROB_TT1[.,5]~(1-p)*PROB_TT1[.,6]~q*PROB_TT1[.,7]~(1-p)*PROB_TT1[.,8]~ (1-q)*PROB_TT1[.,1]~p*PROB_TT1[.,2]~(1-q)*PROB_TT1[.,3]~p*PROB_TT1[.,4]~ (1-q)*PROB_TT1[.,5]~p*PROB_TT1[.,6]~(1-q)*PROB_TT1[.,7]~p*PROB_TT1[.,8]; @CALCULATE MORE PROBABILITIES GIVEN INFORMATION UP TO T@ @Pr[S_t,S_t-1,S_t-k+1,S_t-k,S_t+1*|Y_t]@ @0000 0010 0100 0110 1000 1010 1100 1110 0001 0011 0101 0111 1001 1011 1101 1111 by k-1@ PROB_TT3=zeros(k-1,16); PROB_TT3[.,1:2]=PROB_TT2[.,1:2]; PROB_TT3[2:k-1,3:4]=PROB_TT2[1:k-2,3:4]; PROB_TT3[.,5:6]=PROB_TT2[.,5:6]; PROB_TT3[2:k-1,7:8]=PROB_TT2[1:k-2,7:8]; PROB_TT3[1:k-2,9:10]=PROB_TT2[2:k-1,9:10]; PROB_TT3[.,11:12]=PROB_TT2[.,11:12]; PROB_TT3[1:k-2,13:14]=PROB_TT2[2:k-1,13:14]; PROB_TT3[.,15:16]=PROB_TT2[.,15:16]; @COLLAPSE PROBABILITIES FOR NEXT ITERATION OF HAMILTON FILTER@ @Pr[S_t,S_t-k+1,S_t+1*|Y_t]@ @00 01 10 11 by k-1@ PROB_LL=PROB_TT3[.,1]+PROB_TT3[.,2]+PROB_TT3[.,3]+PROB_TT3[.,4]~ PROB_TT3[.,5]+PROB_TT3[.,6]+PROB_TT3[.,7]+PROB_TT3[.,8]~ PROB_TT3[.,9]+PROB_TT3[.,10]+PROB_TT3[.,11]+PROB_TT3[.,12]~ PROB_TT3[.,13]+PROB_TT3[.,14]+PROB_TT3[.,15]+PROB_TT3[.,16]; likv=likv+ln(f_y); j_iter=j_iter+1; endo; RETP(-LIKV); ENDP; @========================================================================@ PROC (2)=FILTER(PRMTR1); LOCAL PRMTR, q,p,SIG0,SIG1,MU0,MU1,STP,PROB_LL,PROB_TL,s_index, PROB_TT1,PROB_TT2,PROB_TT3,f_ys,likv,j_iter,yt,xt,f_y, beta,ggg,PROBS_TL,PROBS_TT,ddd,MU1STAR; EXTERNAL PROC TRANS, TRANSA, V_PROB; PRMTR=TRANS(PRMTR1|PRMTR_SI'); @PRMTR=TRANSA(PRMTR1);@ @PRMTR=TRANS(PRMTR1);@ SIG0=PRMTR[1]; SIG1=PRMTR[2]; MU0=PRMTR[3]*ones(k-1,1); MU1STAR=PRMTR[4]*ones(k-1,1); DDD=PRMTR[5]; MU1=DDD; GGG=MU1STAR-DDD; beta=PRMTR[6]; q=PRMTR[7]; p=PRMTR[8]; @STP=(1-SUMC(PRMTR[8+1:8+(k-1)*4-1]))|PRMTR[8+1:8+(k-1)*4-1];@ STP=ones(4*(k-1),1)/(4*(k-1)); /* HAMILTON FILTER */ /*==============================================================*/ s_index=cumsumc(ones(k-1,1))-ones(k-1,1); @S_t*=S_t-2 + S_t-3 + ... + S_t-k+1@ @Startup Probability for Hamilton filter@ @Pr[S_t-1,S_t-k,S_t*|Y_t-1]@ @ 00 01 10 11 by k-1 @ PROB_LL=reshape(stp,k-1,4); PROBS_TL=(1-q)/(2-p-q)*ones(T,1); PROBS_TT=zeros(S,1); j_iter=1; do until j_iter>S; yt=yy[j_iter,.]*ones(k-1,1); xt=xx[j_iter,.]*ones(k-1,1); @CALCULATE PROBABILITIES GIVEN INFORMATION UP TO T-1@ @Pr[S_t,S_t-1,S_t-k,S_t*|Y_t-1]=Pr[S_t|S_t-1]*Pr[S_t-1,S_t-k,S_t*|Y_t-1]@ @000 001 010 011 100 101 110 111 by k-1@ PROB_TL= q*PROB_LL[.,1]~q*PROB_LL[.,2]~(1-p)*PROB_LL[.,3]~(1-p)*PROB_LL[.,4]~ (1-q)*PROB_LL[.,1]~(1-q)*PROB_LL[.,2]~p*PROB_LL[.,3]~p*PROB_LL[.,4]; @Pr[S_t=1|Y_t-1]@ PROBS_TL[k+j_iter,1]=sumc(vecr(PROB_TL[.,5:8])); @CALCULATE MARGINAL DENSITY CONDITIONAL ON STATES@ @f(y_t|S_t,S_t-1,S_t-k,S_t*,Y_t-1)@ f_ys=V_PROB(SIG0^2,yt-MU0-PROBS_TL[k+j_iter,1]*GGG- beta*(xt-k*MU0-s_index[.,1].*MU1-sumc(PROBS_TL[j_iter:k-1+j_iter,1])*GGG))~ V_PROB(SIG0^2,yt-MU0-PROBS_TL[k+j_iter,1]*GGG-beta* (xt-k*MU0-(s_index[.,1]+ones(k-1,1)).*MU1- sumc(PROBS_TL[j_iter:k-1+j_iter,1])*GGG))~ V_PROB(SIG0^2,yt-MU0-PROBS_TL[k+j_iter,1]*GGG-beta* (xt-k*MU0-(s_index[.,1]+ones(k-1,1)).*MU1- sumc(PROBS_TL[j_iter:k-1+j_iter,1])*GGG))~ V_PROB(SIG0^2,yt-MU0-PROBS_TL[k+j_iter,1]*GGG-beta* (xt-k*MU0-(s_index[.,1]+2*ones(k-1,1)).*MU1- sumc(PROBS_TL[j_iter:k-1+j_iter,1])*GGG))~ V_PROB(SIG1^2,yt-MU0-MU1-PROBS_TL[k+j_iter,1]*GGG-beta* (xt-k*MU0-s_index[.,1].*MU1-sumc(PROBS_TL[j_iter:k-1+j_iter,1])*GGG))~ V_PROB(SIG1^2,yt-MU0-MU1-PROBS_TL[k+j_iter,1]*GGG-beta* (xt-k*MU0-(s_index[.,1]+ones(k-1,1)).*MU1- sumc(PROBS_TL[j_iter:k-1+j_iter,1])*GGG))~ V_PROB(SIG1^2,yt-MU0-MU1-PROBS_TL[k+j_iter,1]*GGG-beta* (xt-k*MU0-(s_index[.,1]+ones(k-1,1)).*MU1- sumc(PROBS_TL[j_iter:k-1+j_iter,1])*GGG))~ V_PROB(SIG1^2,yt-MU0-MU1-PROBS_TL[k+j_iter,1]*GGG-beta* (xt-k*MU0-(s_index[.,1]+2*ones(k-1,1)).*MU1- sumc(PROBS_TL[j_iter:k-1+j_iter,1])*GGG)); @CALCULATE MARGINAL DENSITY BY SUMMING ACROSS STATES@ @f(y_t|Y_t-1)@ f_y=sumc(vecr(f_ys.*PROB_TL)); @CALCULATE PROBABILITIES GIVEN INFORMATION UP TO T@ @Pr[S_t,S_t-1,S_t-k,S_t*|Y_t]@ @000 001 010 011 100 101 110 111 by k-1@ PROB_TT1=(f_ys.*PROB_TL)/f_y; @CALCULATE MORE PROBABILITIES GIVEN INFORMATION UP TO T@ @Pr[S_t,S_t-1,S_t-k,S_t-k+1,S_t*|Y_t]= Pr[S_t-k+1|S_t-k]*Pr[S_t,S_t-1,S_t-k,S_t*|Y_t]@ @0000 0010 0100 0110 1000 1010 1100 1110 0001 0011 0101 0111 1001 1011 1101 1111 by k-1@ PROB_TT2= q*PROB_TT1[.,1]~(1-p)*PROB_TT1[.,2]~q*PROB_TT1[.,3]~(1-p)*PROB_TT1[.,4]~ q*PROB_TT1[.,5]~(1-p)*PROB_TT1[.,6]~q*PROB_TT1[.,7]~(1-p)*PROB_TT1[.,8]~ (1-q)*PROB_TT1[.,1]~p*PROB_TT1[.,2]~(1-q)*PROB_TT1[.,3]~p*PROB_TT1[.,4]~ (1-q)*PROB_TT1[.,5]~p*PROB_TT1[.,6]~(1-q)*PROB_TT1[.,7]~p*PROB_TT1[.,8]; @CALCULATE MORE PROBABILITIES GIVEN INFORMATION UP TO T@ @Pr[S_t,S_t-1,S_t-k+1,S_t-k,S_t+1*|Y_t]@ @0000 0010 0100 0110 1000 1010 1100 1110 0001 0011 0101 0111 1001 1011 1101 1111 by k-1@ PROB_TT3=zeros(k-1,16); PROB_TT3[.,1:2]=PROB_TT2[.,1:2]; PROB_TT3[2:k-1,3:4]=PROB_TT2[1:k-2,3:4]; PROB_TT3[.,5:6]=PROB_TT2[.,5:6]; PROB_TT3[2:k-1,7:8]=PROB_TT2[1:k-2,7:8]; PROB_TT3[1:k-2,9:10]=PROB_TT2[2:k-1,9:10]; PROB_TT3[.,11:12]=PROB_TT2[.,11:12]; PROB_TT3[1:k-2,13:14]=PROB_TT2[2:k-1,13:14]; PROB_TT3[.,15:16]=PROB_TT2[.,15:16]; @COLLAPSE PROBABILITIES FOR NEXT ITERATION OF HAMILTON FILTER@ @Pr[S_t,S_t-k+1,S_t+1*|Y_t]@ @00 01 10 11 by k-1@ PROB_LL=PROB_TT3[.,1]+PROB_TT3[.,2]+PROB_TT3[.,3]+PROB_TT3[.,4]~ PROB_TT3[.,5]+PROB_TT3[.,6]+PROB_TT3[.,7]+PROB_TT3[.,8]~ PROB_TT3[.,9]+PROB_TT3[.,10]+PROB_TT3[.,11]+PROB_TT3[.,12]~ PROB_TT3[.,13]+PROB_TT3[.,14]+PROB_TT3[.,15]+PROB_TT3[.,16]; @COLLAPSE PROBABILITIES TO GET FILTERED PROB OF S_t|t@ @Pr[S_t=1|Y_t]@ PROBS_TT[j_iter,1]=sumc(vecr(PROB_LL[.,3:4])); j_iter=j_iter+1; endo; RETP(PROBS_TL[k+1:T],PROBS_TT); ENDP; @>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@ PROC TRANS(coef0); @ constraining values of reg. coeff.@ local coef1; coef1=coef0; coef1[1:2]=exp(-coef0[1:2]); coef1[3]=exp(-coef0[3])/(1+exp(-coef0[3])); coef1[4:5]=0.01*coef0[4:5]; coef1[7:8]=exp(-coef0[7:8])./(1+exp(-coef0[7:8])); coef1[8+1:8+4*(k-1)-1,1]= ((coef0[8+1:8+4*(k-1)-1,1])^2) /(1+sumc(coef0[8+1:8+4*(k-1)-1,1].^2)); retp(coef1); endp; @==========================================================@ @>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@ PROC TRANSA(coef0); @ constraining values of reg. coeff.@ local coef1; coef1=coef0; coef1[1:2]=exp(-coef0[1:2]); coef1[3]=exp(-coef0[3])/(1+exp(-coef0[3])); coef1[4:5]=0.01*coef0[4:5]; coef1[7:8]=exp(-coef0[7:8])./(1+exp(-coef0[7:8])); retp(coef1); endp; @==========================================================@ @>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@ PROC V_PROB(HE, EV); @ CALCULATES Pr[Yt/St,Yt-1]@ LOCAL VAL; VAL=(1/SQRT(2*PI*HE))*EXP(-0.5*EV.*EV/HE); RETP(VAL); ENDP; @>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@ PROC SMOOTH(pr_tl1,pr_tt1); @pr_TT1 contains Pr[S_t|Y_t]@ @pr_TL1 contains Pr[S_t|Y_{t-1}]@ local ppr, qpr, pr_sm0,pr_sm1, j_iter,pr_sm00,pr_sm01,pr_sm10,pr_sm11, pr_tt0,pr_tl0,TT; TT=ROWS(PR_TT1); PPR=PRM_fnl[8]; @Pr[St=1/St-1=1]@ QPR=PRM_fnl[7]; @Pr[St=0/St-1=0]@ pr_tt0=1-pr_tt1; pr_tl0=1-pr_tl1; pr_sm0=pr_tt0; @ pr_sm0 will contain Pr[S_t|Y_T]@ pr_sm1=pr_tt1; j_iter=TT-1; do until j_iter < 1; @The followings are P[S_t, S_t+1|Y_T] @ pr_sm00=pr_sm0[j_iter+1,1]*qpr* pr_tt0[j_iter,1]/ pr_tl0[j_iter+1,1]; pr_sm01=pr_sm1[j_iter+1,1]*(1-qpr)*pr_tt0[j_iter,1]/ pr_tl1[j_iter+1,1]; pr_sm10=pr_sm0[j_iter+1,1]*(1-ppr)*pr_tt1[j_iter,1]/ pr_tl0[j_iter+1,1]; pr_sm11=pr_sm1[j_iter+1,1]*ppr* pr_tt1[j_iter,1]/ pr_tl1[j_iter+1,1]; pr_sm0[j_iter,1]=pr_sm00+pr_sm01; pr_sm1[j_iter,1]=pr_sm10+pr_sm11; j_iter=j_iter -1; endo; RETP(pr_sm1); @This proc returns Pr[S_t=1|Y_T]@ endp; @++++++++++++++++++++++++++++++++++++++++++++++++++++++++@