/* April. 6, 2009 This program is designed to calculate the model-free measure of the business cycle for US real GDP For more details, see "The Asymmetric Business Cycle" and "Practical Computation of the Modle-Free Business Cycle", available at http://artsci.wustl.edu/~morley/ To use this program, save the data files gdp4706.txt, gdp4708.txt, and any other data file of interest into the working directory with this file. Also, create a folder "output" in the working directory. The various estimates and data output files will be stored there. Note: the output file lik.out stores the following info sample size~model-type~AR order~Error dist~sb assumption~loglik~r# of parameters~AIC~BIC~MSC; where for model type: AR=0, BBU=1, BBV=2, BBD=3 sb assumption: no break=1, scaled break=8 The weights for the model free measure are AR(2) 25% AR(12) 25% BBD-0 50% The mean of the model-free measure is estimated using the BB model (the other models assume mean=0) */ new; cls; library optmum,pgraph; format /m1 /rd 7,2; @GLOBALS for ALL MODELS@ error_dist=1; @1=normal, 2=t-dist@ sb=8; @1=no break, 8=mb+vb as scale shifts@ sample=1; @set=1 if updated sample, 0 for original sample@ @DATA@ if sample==0; load data[240,2]=gdp4706.txt; @Real GDP 1947:Q1 -- 2006:Q4; @ endif; if sample==1; load data[248,2]=gdp4708.txt; @Real GDP 1947:Q1 -- 2008:Q4; (final release)@ endif; @OUTPUT FILES@ output file=output/lik.out reset; output file=output/prmtr.out reset; output file=output/bn2.dta reset; output file=output/bn12.dta reset; output file=output/bbd0.dta reset; output file =output/probabilities.txt reset; output file=output/model_free.dta reset; output off; @AR MODELS ESTIMATION@ @Globals for AR models to loop over@ p=2; @AR order-- program jumps from p=2 to p=12@ m=0; @Model type: m=0 => AR@ y_level_t = ln(data[.,2])*100; maxp=12; @maximum number of AR lags; code set up for anything up to 12@ mu_=(y_level_t[rows(y_level_t)]-y_level_t[1])/(rows(y_level_t)-1); j=1; do until j>maxp; y_level_t=y_level_t[1]-mu_|y_level_t; j=j+1; endo; sdate=maxp+2; @1947:1 for differences@ yy = y_level_t[sdate:rows(y_level_t)] - y_level_t[sdate-1:rows(y_level_t)-1]; x_mat_full={}; j=1; do until j>maxp; x_mat_full = x_mat_full~(y_level_t[sdate-j:rows(y_level_t)-j] - y_level_t[sdate-(j+1):rows(y_level_t)-(j+1)]); j=j+1; endo; t = rows(yy); @Structural Break dates@ base_year=1947; year1=1973; quarter1=1; year2=1984; quarter2=2; breakdate1=4*(year1-base_year)+quarter1-(sdate-maxp-1); @1973:Q1, meaning new parameter in Q2@ breakdate2=4*(year2-base_year)+quarter2-(sdate-maxp-1); @1984:Q2, meaning new parameter in Q3@ output file=output/prmtr.out; output on; "AR Models (p~error_dist~sb)"; "Estimated parameters"; "(standard errors)"; ""; output off; counter=0; count_lik=0; do until p>12; x_mat=x_mat_full[.,1:p]; prmtr_in=zeros(2+p,1); prmtr_in[2]=meanc(yy); b=yy/(ones(rows(x_mat),1)~x_mat); prmtr_in[3:2+p]=b[2:1+p]; prmtr_in[1]=ln(stdc(yy-(ones(rows(x_mat),1)~x_mat)*b)^2); if error_dist==2; prmtr_in=prmtr_in|-10; endif; if sb==8; prmtr_in=prmtr_in|0|0; endif; @ Maximum Likelihood Estimation @ @==================================================@ {xout,fout,gout,cout}=optmum(&lik_fcn_ar,PRMTR_in); PRM_FNL=TRANS_ar(xout); "Calculating Hessian"; hout0=hessp(&lik_fcn_ar,xout); hout=inv(hout0); grdn_fnl=gradfd(&TRANS_ar,xout); Hsn_fnl=grdn_fnl*hout*grdn_fnl'; SD_fnl =sqrt(diag(Hsn_fnl)); t_1=t; k_s=0; n__=1; @number of regimes@ lambda_1=n__; lambda_2=n__; delta_1=1; delta_2=1; msc_penalty=-0.5*t + (rows(prm_fnl)-n__*(n__-1)-n__*k_s-1)*t / maxc(( t-(rows(prm_fnl)-n__*(n__-1)-n__*k_s-1)-1)|1) + 0.5*t_1*(t_1+lambda_1*k_s)/maxc((delta_1*t_1 - lambda_1*k_s -2)|1); output file = output/lik.out; output on; t~m~p~error_dist~sb~-fout~rows(prm_fnl)~(-fout - rows(prm_fnl))~(-fout - rows(prm_fnl)*0.5*ln(t))~(-fout-msc_penalty); output off; output file = output/prmtr.out; output on; p~error_dist~sb; prm_fnl'; "(";;sd_fnl';;")"; ""; output off; @BN Decomposition@ mu_mat=prm_fnl[2]*ones(p+rows(yy),1); @unconditional mean@ if sb==8; if error_dist==1; mu_new=prm_fnl[2+p+1]*prm_fnl[2]; else; mu_new=prm_fnl[2+p+2]*prm_fnl[2]; endif; mu_mat[p+breakdate1+1:p+rows(yy)]=mu_new*ones(rows(yy)-breakdate1,1); endif; if p==1; B_t=yy-mu_mat[p+1:rows(mu_mat)]; else; B_t=yy-mu_mat[p+1:rows(mu_mat)]; i=1; do until i>p-1; B_t=B_t~x_mat[.,i]-mu_mat[p+1-i:rows(mu_mat)-i]; i=i+1; endo; endif; if p==1; F=prm_fnl[3]; else; F= prm_fnl[3:2+p,1]'| eye(p-1)~zeros(p-1,1); @companion matrix for AR(p)@ endif; c_tt=(-F*inv(eye(p)-F)*B_t')'; @BN cycle@ c_mat={}; if p==2; output file=output/bn2.dta; output on; c_tt[.,1]; output off; bn2=c_tt[.,1]; endif; if p==12; output file=output/bn12.dta; output on; c_tt[.,1]; output off; bn12=c_tt[.,1]; endif; p=p+10; endo; @BOUNCEBACK MODEL ESTIMATION@ @Globals for BB model@ m=3; @model type (1=BBU, 2=BBV, 3=BBD)@ p=0; @ar lags@ norm_state=0; @NORMALIZE WHICH REGIME IS ZERO-MEAN STATE@ y_level_t = ln(data[.,2])*100; maxp=2; @maximum number of AR lags; code set up for anything up to 12@ MSTAR = 6; mu_=(y_level_t[rows(y_level_t)]-y_level_t[1])/(rows(y_level_t)-1); maxp=2; j=1; do until j>maxp+mstar; y_level_t=y_level_t[1]-mu_|y_level_t; j=j+1; endo; sdate=maxp+mstar+2; y_level=y_level_t[sdate:rows(y_level_t)]; yy = y_level_t[sdate:rows(y_level_t)] - y_level_t[sdate-1:rows(y_level_t)-1]; x_mat_full={}; j=1; do until j>maxp+mstar; x_mat_full = x_mat_full~(y_level_t[sdate-j:rows(y_level_t)-j] - y_level_t[sdate-(j+1):rows(y_level_t)-(j+1)]); j=j+1; endo; t = rows(yy); @Structural Break dates@ base_year=1947; year1=1973; quarter1=1; year2=1984; quarter2=2; breakdate1=4*(year1-base_year)+quarter1-(sdate-maxp-mstar-1); @1973:Q1, meaning new parameter in Q2@ breakdate2=4*(year2-base_year)+quarter2-(sdate-maxp-mstar-1); @1984:Q2, meaning new parameter in Q3@ output file=output/prmtr.out; output on; "Bounceback Model (p~error_dist~sb)"; "Estimated parameters"; "(standard errors)"; ""; output off; count_lik=0; ctt_mat={}; counter=0; NO_ST=p+mstar+1; @ NUMBER OF STATES TO BE CONSIDERED@ DMNSION = 2^NO_ST; @ NUMBER OF CASES TO BE CONSIDERED@ st_mat=zeros(DMNSION,NO_ST); j=1; if p==0; st6=0;do until st6>1; st5=0;do until st5>1; st4=0;do until st4>1; st3=0;do until st3>1; st2=0;do until st2>1; st1=0;do until st1>1; st=0;do until st>1; st_mat[j,.] = st6~st5~st4~st3~st2~st1~st; j=j+1; st=st+1;endo; st1=st1+1;endo; st2=st2+1;endo; st3=st3+1;endo; st4=st4+1;endo; st5=st5+1;endo; st6=st6+1;endo; endif; x_mat=x_mat_full[.,1:p+mstar]; @STARTING VALUES for ESTIMATION@ n1=5; @initial # of different sets of starting values@ _opmiter=3*(1e+2); @max # of iterations for optimization@ prmtr_in1=zeros(n1,6+p); prm_fnl_vec1={}; xout_vec={}; prmtr_in1n={}; jjj=1; do until jjj>n1; prmtr_in1[jjj,1]=3.2+rndn(1,1); prmtr_in1[jjj,2]=1.2+rndn(1,1); prmtr_in1[jjj,3]=0.8+0.2*rndn(1,1); prmtr_in1[jjj,4]=0.8; prmtr_in1[jjj,5]=rndn(1,1); prmtr_in1[jjj,6]=-0.3+0.5*rndn(1,1); PRMTR_IN=prmtr_in1[jjj,.]'; if sb==8; prmtr_in=prmtr_in|0|0; endif; if error_dist==2; prmtr_in=prmtr_in|-10; endif; prmtr_in1n=prmtr_in1n|prmtr_in'; /* Maximum Likelihood Estimation */ {xout,fout,gout,cout} = optmum(&lik_fcn,PRMTR_in); prm_fnl_vec1=prm_fnl_vec1|(trans(xout)'~-fout); xout_vec=xout_vec|(xout'~-fout); jjj=jjj+1; endo; prmtr_in=prmtr_in1n~xout_vec[.,cols(xout_vec)]; prmtr_in=rev(sortc(prmtr_in,rows(xout)+1)); prm_fnl_vec1=rev(sortc(prm_fnl_vec1,rows(xout)+1)); prm_fnl=prm_fnl_vec1[1,1:rows(xout)]'; fout=-prm_fnl_vec1[1,rows(xout)+1]; xout_vec=rev(sortc(xout_vec,rows(xout)+1)); xout=xout_vec[1,1:rows(xout)]'; "Calculating Hessian..... Please be patient!!!!"; hout0=hessp(&lik_fcn,xout); hout=inv(hout0); grdn_fnl=gradp(&TRANS,xout); Hsn_fnl=grdn_fnl*hout*grdn_fnl'; SD_fnl =sqrt(diag(Hsn_fnl)); @Standard errors of the estimated coefficients@ /* Calculate and Print Regime Probabilities */ {pr_tt0,pr_tl0} = FILTER(XOUT); @Pr[S_t=0|Y_t] and Pr[S_t=0|Y_{t-1}]@ smooth0 = smooth(xout); @Pr[S_t=0|Y_T], Smoothed probabilities@ output file = output/probabilities.txt; output on; (1-pr_tt0)~(1-smooth0); output off; t_1=sumc(smooth0); t_2=sumc(1-smooth0); lambda_1=2; lambda_2=2; delta_1=1; delta_2=1; k_s=1; n__=2; @number of regimes@ msc_penalty=-0.5*t + (rows(prm_fnl)-n__*(n__-1)-n__*k_s-1)*t / maxc(( t-(rows(prm_fnl)-n__*(n__-1)-n__*k_s-1)-1)|1) + 0.5*t_1*(t_1+lambda_1*k_s)/maxc((delta_1*t_1 - lambda_1*k_s -2)|1) + 0.5*t_2*(t_2+lambda_2*k_s)/maxc((delta_2*t_2 - lambda_2*k_s -2)|1); output file = output/lik.out; output on; t~m~p~error_dist~sb~-fout~rows(prm_fnl)~(-fout - rows(prm_fnl))~(-fout - rows(prm_fnl)*0.5*ln(t))~(-fout-msc_penalty); output off; output file = output/prmtr.out; output on; p~error_dist~sb; prm_fnl'; "(";;sd_fnl';;")"; ""; output off; @Estimate Trend using RDSS Approach @ p00_fnl=prm_fnl[1]; @Pr[St=0/St-1=0]@ p11_fnl=prm_fnl[2]; @Pr[St=1/St-1=1]@ pi_bar_fnl=(1-p00_fnl)/(2-p11_fnl-p00_fnl); @Pr[St=1]@ var_fnl=prm_fnl[3]^2; mu_bar_fnl=prm_fnl[4]; d_gam_fnl=prm_fnl[5]; lambda_fnl=prm_fnl[6]; gam0_fnl=mu_bar_fnl; @ recession, vs. boom@ gam1_fnl=d_gam_fnl; @ recession, vs. boom@ if p==0; phi_fnl=0; else; phi_fnl= prm_fnl[7:6+p]; endif; @Structural break parameters@ d_gam_pb_fnl=d_gam_fnl; mu_bar_pb_fnl=mu_bar_fnl; mu_bar_pb2_fnl=mu_bar_fnl; var_pb_fnl=var_fnl; if sb==8; mu_bar_pb_fnl=prm_fnl[6+p+1]*mu_bar_fnl; mu_bar_pb2_fnl=mu_bar_pb_fnl; d_gam_pb_fnl=prm_fnl[6+p+2]*d_gam_fnl; var_pb_fnl=prm_fnl[6+p+2]^2*var_fnl; endif; trend_est = {}; if p==0; F=0; H=0; elseif p==1; F=phi_fnl[1]; H=1; else; F=phi_fnl[1:p]'|eye(p-1)~zeros(p-1,1); H=1~zeros(1,p-1); endif; qqq=1; do until qqq > T; if qqq>breakdate1; gam0_fnl=mu_bar_pb_fnl; endif; if qqq>breakdate2; gam1_fnl=d_gam_pb_fnl; gam0_fnl=mu_bar_pb2_fnl; endif; lr_vec = {}; p_iter = 1; do until p_iter > dmnsion; s_vec = st_mat[p_iter,.]; s_vec = s_vec'; @fill with S_t-m-p ... S_t-1 S_t@ ydiff=zeros(mstar+p,1); @fill with Dy_t ... Dy_t-m-p+1@ ydiff[1] = yy[qqq]; @Dy_t@ j=2; do until j>mstar+p; ydiff[j] = x_mat[qqq,j-1]; j=j+1; endo; ydiff_gamma={}; @fill with (gam1+Dy_t-m+1) ... (gam1+Dy_t)@ j=mstar; do until j<1; if m==1 or m==2; ydiff_gamma=ydiff_gamma|gam1_fnl; elseif m==3; ydiff_gamma=ydiff_gamma|ydiff[j]+gam1_fnl; endif; j=j-1; endo; if p ge 1; y_tilda={}; @fill with Dy_t ... Dy_t-p+1@ mu_tilda={}; @fill with mu_t ... mu_t-p+1@ j=1; do until j>p; y_tilda=y_tilda|ydiff[j]; if m==1 or m==3; mu_tilda=mu_tilda| gam0_fnl + gam1_fnl*s_vec[mstar+p+1+1-j] + lambda_fnl*(s_vec[p+1+1-j:p+mstar+1-j]'*ydiff_gamma); elseif m==2; mu_tilda=mu_tilda| gam0_fnl + gam1_fnl*s_vec[mstar+p+1+1-j] + lambda_fnl*(1-s_vec[mstar+p+1+1-j])*(s_vec[p+1+1-j:p+mstar+1-j]'*ydiff_gamma); endif; j=j+1; endo; endif; nl_momentum=0; f_iter=1; do until f_iter>mstar; nl_momentum=nl_momentum+(mstar-f_iter+1)*lambda_fnl*ydiff_gamma[mstar-f_iter+1]*s_vec[p+1+mstar-f_iter+1]; @m*(gam1+Dy_t)*S_t + (m-1)*(gam1+Dy_t-1)*S_t-1 + ...@ f_iter=f_iter+1; endo; if p==0; lr_fcast = y_level[qqq] + nl_momentum; else; lr_fcast = y_level[qqq] + H*F*inv(eye(rows(F))-F)*(y_tilda-mu_tilda) + nl_momentum; endif; lr_vec = lr_vec|lr_fcast; p_iter = p_iter+1; endo; prob_mat = get_probs(xout,qqq); trend_est = trend_est|sumc(lr_vec.*prob_mat); qqq=qqq+1; endo; @RDSS Estimate of Cycle@ c_tt = y_level-trend_est; ctt_mat=ctt_mat~c_tt; output file=output/bbd0.dta; output on; c_tt; output off; bbd0=c_tt; @CONSTRUCT MODEL-FREE CYCLE@ weights_ ={ 0.25 0.25 0.5 }; weights_=weights_'; rdss_all= bn2~ bn12~ bbd0; weights_=weights_./sumc(weights_); rdss_free = sumc((weights_.*rdss_all')); rdss_free=(meanc(bbd0[.,1])/meanc(rdss_free))*rdss_free; output file=output/model_free.dta; output on; rdss_free~bbd0; output off; begwind; window(3,1,0); setwind(1); title("Model-Free Cycle"); xy(rows(rdss_free),rdss_free~zeros(rows(rdss_free),1)); nextwind; title("BBD0 Cycle"); xy(rows(rdss_free),bbd0~zeros(rows(rdss_free),1)); nextwind; title("Smoothed Probability of Recession State for BBD0 Model"); xy(rows(rdss_free),(1-smooth0)); endwind; cls; load comparison[3,10] = output/lik.out; "sample size~model-type~AR order~Error dist~sb assumption~loglik~r# of parameters~AIC~BIC~MSC" comparison; ""; "where for model type: AR=0, BBU=1, BBV=2, BBD=3 and for sb assumption: no break=1, scaled break=8"; ""; "Weights for Model-Free Measure" "AR(2)~AR(12)~BBD"; weights_'; end; @========================================================================@ proc LIK_FCN_ar(PRMTR1); local prmtr, lik, j_iter, f_cast1, pr_vl, likv, phi, j, MU, VAR,NU_HAT,var_pb,mu_pb; PRMTR=TRANS_ar(PRMTR1); VAR = PRMTR[1]; MU = PRMTR[2]; PHI = prmtr[3:2+p]; if error_dist==2; if PRMTR[2+p+1]<1000; NU_HAT=PRMTR[2+p+1]; else; nu_hat=1000; endif; endif; var_pb=var; mu_pb=mu; if sb==8; if error_dist==1; mu_pb=prmtr[2+p+1]*mu; var_pb=prmtr[2+p+2]^2*var; else; mu_pb=prmtr[2+p+2]*mu; var_pb=prmtr[2+p+3]^2*var; endif; endif; LIKV=0.0; J_ITER=1; do until J_ITER>T; if j_iter>breakdate1; mu=mu_pb; endif; if j_iter>breakdate2; var=var_pb; endif; F_CAST1 = YY[J_ITER,1] -mu -(X_MAT[J_ITER,.]-mu)*PHI; if error_dist==1; PR_VL=(1./SQRT(2.*PI.*VAR)).*EXP(-0.5*F_CAST1.*F_CAST1./VAR); else; PR_VL=studt_d(f_cast1,0,VAR,NU_HAT); endif; LIK=LN(PR_VL); LIKV = LIKV+LIK; J_ITER = J_ITER+1; endo; retp(-LIKV); endp; @========================================================================@ @================================================================@ proc TRANS_ar(c0); @ constraining values of reg. coeff.@ local c1; c1 = c0; c1[1]=exp(c0[1]); if error_dist==2; c1[2+p+1]=1+10000*(exp(c0[2+p+1])/(1+exp(c0[2+p+1]))); endif; if sb==8; if error_dist==1; c1[2+p+1]=exp(c0[2+p+1]); c1[2+p+2]=exp(c0[2+p+2]); else; c1[2+p+2]=exp(c0[2+p+2]); c1[2+p+3]=exp(c0[2+p+3]); endif; endif; retp(c1); endp; @========================================================================@ proc STUDT_D(MUFNL0,MU0,V0,NU0); local VAL,lnnu,lnnup; if (nu0/2)<=1; lnnu=ln(gamma(nu0/2)); else; lnnu=lnfact((nu0/2)-1); endif; if ((nu0+1)/2)<=1; lnnup=ln(gamma((nu0+1)/2)); else; lnnup=lnfact(((nu0+1)/2)-1); endif; @Scalar version@ VAL= lnnup - (1/2)*ln(v0*nu0*pi) - lnnu - ((nu0+1)/2)*ln(1+(1/nu0)*(MUFNL0-MU0).*(MUFNL0-MU0)./v0); retp(exp(VAL)); endp; @==================================================================@ @========================================================================@ proc LIK_FCN(PRMTR1); local prmtr, ppr,qpr,prob__0,prob__1,QQ, lik, j_iter, pr__0_l,pr__1_l, var_L, pr_val,likv,phi,PSIX, vecp,st,st1,st2,st3,st4, ST5,ST6,ST7,ST8,ST9,ST10,ST11,ST12,ST13,pr_tr,pr_trf1,pr_trf,prob__t, prob__,pro_,pr_vl,j,psi1,psi2,var_c,delta0,DELTA1,MU0,MU1,st_k,st_k1,st_k2, st_k3,st_k4,f_cast1,f_cast2,PR_VL1,PR_VL2,pr_trf0,PR_TRF2,PR_TRF3,PR_TRF4, PR_TRF5,PR_TRF6,PR_TRF7,PR_TRF8,psic,psiL,TMPRY1,TMPRY2,SM_PRL,TMP_P0, JJJ,MU_MAT,D_MAT,FLT_PRN,F1,F2,TMP00,TMP0,SM_PR0,SM_PR00,prob_dd,VAR,A,EN,omga1, omga2,omga3,omga4,cdr1,cdr2,cdr3,cdr4,cdr1_m,cdr2_m,cdr3_m,cdr4_m,var1,var2, mu01,mu11,mu02,mu12,mu_add,pi_bar,mu_bar,d_mu,mu_pb,d_mu_pb,p00,p11,d_gam,lambda,gam0,gam1, mu_bar_pb,var_pb,d_gam_pb,i_iter,nu_hat,mu_bar_pb2; PRMTR=TRANS(PRMTR1); p00=prmtr[1]; @Pr[St=0/St-1=0]@ p11=prmtr[2]; @Pr[St=1/St-1=1]@ pi_bar=(1-p00)/(2-p11-p00); @Pr[St=1]@ var=prmtr[3]^2; mu_bar=prmtr[4]; d_gam=prmtr[5]; lambda=prmtr[6]; gam0=mu_bar; gam1=d_gam; @ recession, vs. boom@ if p==0; phi=0; else; phi= prmtr[7:6+p]; endif; @Structural break parameters@ d_gam_pb=d_gam; mu_bar_pb=mu_bar; mu_bar_pb2=mu_bar; var_pb=var; if sb==8; mu_bar_pb=prmtr[6+p+1]*mu_bar; mu_bar_pb2=mu_bar_pb; d_gam_pb=prmtr[6+p+2]*d_gam; var_pb=prmtr[6+p+2]^2*var; endif; if error_dist==2; if PRMTR[rows(prmtr)]<1000; NU_HAT=PRMTR[rows(prmtr)]; else; nu_hat=1000; endif; endif; MU_MAT = ST_MAT*gam1 + (ONES(DMNSION,NO_ST))*gam0; PR_TR=(p00~ (1-p11))| ((1-p00)~ p11); A = (eye(2)-pr_tr)|ones(1,2); EN=(0|0|1); if det(A'A) le 0; PROB__T=0.5|0.5; else; PROB__T = INV(A'A)*A'EN; @PR[S_t=0]|PR[S_t=1], 2x1 steady-state PROBABILITIES@ endif; PR_TRF0=VEC(PR_TR); i_iter=1; do until i_iter>mstar+p-1; PROB__T=VECR(PROB__T~PROB__T).*PR_TRF0; PR_TRF0=PR_TRF0|PR_TRF0; i_iter=i_iter+1; endo; PR_TRF =PR_TRF0; PROB__=VECR(PROB__T~PROB__T); LIKV=0; J_ITER=1; do until J_ITER>T; if j_iter>breakdate1; gam0=mu_bar_pb; MU_MAT=ST_MAT*gam1 + (ONES(DMNSION,NO_ST))*gam0; endif; if j_iter>breakdate2; VAR=var_pb; gam0=mu_bar_pb2; gam1=d_gam_pb; MU_MAT=ST_MAT*gam1 + (ONES(DMNSION,NO_ST))*gam0; endif; if p==0; if m==1; F_CAST1=(YY[J_ITER,1])*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(st_mat[.,p+1:mstar+p]*rev((gam1*ones(1,mstar))')); elseif m==2; F_CAST1=(YY[J_ITER,1])*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(1-st_mat[.,mstar+p+1]).*(st_mat[.,p+1:mstar+p]*rev((gam1*ones(1,mstar))')); elseif m==3; F_CAST1=(YY[J_ITER,1])*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(st_mat[.,p+1:mstar+p]*rev((x_mat[j_iter,1:mstar]+gam1*ones(1,mstar))')); endif; else; if m==1; F_CAST1=(YY[J_ITER,1]-X_MAT[J_ITER,1:p]*PHI)*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(st_mat[.,p+1:mstar+p]*rev((gam1*ones(1,mstar))')); i_iter=1; do until i_iter>p; F_CAST1=F_CAST1 + MU_MAT[.,p+1-i_iter]*PHI[i_iter] + lambda*(st_mat[.,p+1-i_iter:mstar+p-i_iter]*rev((gam1*ones(1,mstar))'))*PHI[i_iter]; i_iter=i_iter+1; endo; elseif m==2; F_CAST1=(YY[J_ITER,1]-X_MAT[J_ITER,1:p]*PHI)*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(1-st_mat[.,mstar+p+1]).*(st_mat[.,p+1:mstar+p]*rev((gam1*ones(1,mstar))')); i_iter=1; do until i_iter>p; F_CAST1=F_CAST1 + MU_MAT[.,p+1-i_iter]*PHI[i_iter] + lambda*(1-st_mat[.,mstar+1-i_iter]).*(st_mat[.,p+1-i_iter:mstar+p-i_iter]*rev((gam1*ones(1,mstar))'))*PHI[i_iter]; i_iter=i_iter+1; endo; elseif m==3; F_CAST1=(YY[J_ITER,1]-X_MAT[J_ITER,1:p]*PHI)*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(st_mat[.,p+1:mstar+p]*rev((x_mat[j_iter,1:mstar]+gam1*ones(1,mstar))')); i_iter=1; do until i_iter>p; F_CAST1=F_CAST1 + MU_MAT[.,p+1-i_iter]*PHI[i_iter] + lambda*(st_mat[.,p+1-i_iter:mstar+p-i_iter]*rev((x_mat[j_iter,1+i_iter:mstar+i_iter]+gam1*ones(1,mstar))'))*PHI[i_iter]; i_iter=i_iter+1; endo; endif; endif; VAR_L=VAR*ONES(DMNSION,1); PROB_DD=PR_TRF .* PROB__; /* Pr[S_t,S_{t-1},.....,S_{t-5} | I(t-1)]*/ if error_dist==1; PR_VL=(1./SQRT(2.*PI.*VAR)).*EXP(-0.5*F_CAST1.*F_CAST1./VAR).*PROB_DD; else; PR_VL=studt_d(f_cast1,0,VAR,NU_HAT).*PROB_DD; endif; /* 2^6x1 */ /* Joint density of y_t,S_t,..,S_{t-5} given past information : */ PR_VAL=SUMC(PR_VL); /* f(y_t|I(t-1)), density of y_t given past information: This is weighted average of 2^6 conditional densities */ LIK=-1*LN(PR_VAL); PRO_=PR_VL/PR_VAL; /* Pr[S_t,S_{t-1},...,S_{t-5} | I(t-1),y_t]*/ /* Updating of prob. with new information, y_t */ PROB__T=PRO_[1:DMNSION/2,1]+PRO_[DMNSION/2+1:DMNSION,1]; /* Integrate out S_{t-4}: then you get Pr[S_t, S_{t-1},....,S_{t-4}| Y_t] */ /* 2^5x1*/ PROB__=VECR(PROB__T~PROB__T); /* 2^6x1 */ LIKV = LIKV+LIK; J_ITER = J_ITER+1; endo; retp(LIKV); endp; /**********************************************************************/ proc (2) = FILTER(PRMTR1); local prmtr, ppr,qpr,prob__0,prob__1,QQ, lik, j_iter, pr__0_l,pr__1_l, F_cast, var_L, pr_val,likv,phi,PSIX, vecp,st,st1,st2,st3,st4,ST5,ST6,ST7,ST8,ST9,ST10,ST11,ST12,ST13, pr_tr,pr_trf1,pr_trf,prob__t,prob__,pro_,pr_vl,j,psi1,psi2,var_c, delta0,DELTA1,MU0,MU1,st_k,st_k1,st_k2,st_k3,st_k4, f_cast1,f_cast2,PR_VL1,PR_VL2,pr_trf0, PR_TRF2,PR_TRF3,PR_TRF4,PR_TRF5,PR_TRF6,PR_TRF7,PR_TRF8,psic,psiL,PR_STL0,PR_STT0, TMPRY1,TMPRY2,SM_PRL,TMP_P0,JJJ,MU_MAT,D_MAT,FLT_PRN, F1,F2,TMP00,TMP0,SM_PR0,SM_PR00,prob_dd,VAR,TMP,OUT_MAT,A,EN,omga1,omga2, omga3,omga4,cdr1,cdr2,cdr3,cdr4,cdr1_m,cdr2_m,cdr3_m,cdr4_m,var1,var2, mu01,mu11,mu02,mu12,mu_add,pi_bar,mu_bar,d_mu,mu_pb,d_mu_pb,p00,p11,d_gam,lambda,gam0,gam1, mu_bar_pb,var_pb,d_gam_pb,i_iter,nu_hat,mu_bar_pb2; PRMTR=TRANS(PRMTR1); p00=prmtr[1]; @Pr[St=0/St-1=0]@ p11=prmtr[2]; @Pr[St=1/St-1=1]@ pi_bar=(1-p00)/(2-p11-p00); @Pr[St=1]@ var=prmtr[3]^2; mu_bar=prmtr[4]; d_gam=prmtr[5]; lambda=prmtr[6]; gam0=mu_bar; gam1=d_gam; @ recession, vs. boom@ if p==0; phi=0; else; phi= prmtr[7:6+p]; endif; @Structural break parameters@ d_gam_pb=d_gam; mu_bar_pb=mu_bar; mu_bar_pb2=mu_bar; var_pb=var; if sb==8; mu_bar_pb=prmtr[6+p+1]*mu_bar; mu_bar_pb2=mu_bar_pb; d_gam_pb=prmtr[6+p+2]*d_gam; var_pb=prmtr[6+p+2]^2*var; endif; if error_dist==2; if PRMTR[rows(prmtr)]<1000; NU_HAT=PRMTR[rows(prmtr)]; else; nu_hat=1000; endif; endif; MU_MAT = ST_MAT*gam1 + (ONES(DMNSION,NO_ST))*gam0; PR_TR=(p00~ (1-p11))| ((1-p00)~ p11); A = (eye(2)-pr_tr)|ones(1,2); EN=(0|0|1); if det(A'A) le 0; PROB__T=0.5|0.5; else; PROB__T = INV(A'A)*A'EN; @PR[S_t=0]|PR[S_t=1], 2x1 steady-state PROBABILITIES@ endif; PR_TRF0=VEC(PR_TR); i_iter=1; do until i_iter>mstar+p-1; PROB__T=VECR(PROB__T~PROB__T).*PR_TRF0; PR_TRF0=PR_TRF0|PR_TRF0; i_iter=i_iter+1; endo; PR_TRF =PR_TRF0; PROB__=VECR(PROB__T~PROB__T); PR_STT0=ZEROS(T,1); @WILL save Pr[S_t=0|Y_{t}@ PR_STL0=ZEROS(T,1); @WILL save Pr[S_t=0|Y_{t-1}@ LIKV=0.0; J_ITER=1; do until J_ITER>T; if j_iter>breakdate1; gam0=mu_bar_pb; MU_MAT=ST_MAT*gam1 + (ONES(DMNSION,NO_ST))*gam0; endif; if j_iter>breakdate2; VAR=var_pb; gam0=mu_bar_pb2; gam1=d_gam_pb; MU_MAT=ST_MAT*gam1 + (ONES(DMNSION,NO_ST))*gam0; endif; if p==0; if m==1; F_CAST1=(YY[J_ITER,1])*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(st_mat[.,p+1:mstar+p]*rev((gam1*ones(1,mstar))')); elseif m==2; F_CAST1=(YY[J_ITER,1])*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(1-st_mat[.,mstar+p+1]).*(st_mat[.,p+1:mstar+p]*rev((gam1*ones(1,mstar))')); elseif m==3; F_CAST1=(YY[J_ITER,1])*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(st_mat[.,p+1:mstar+p]*rev((x_mat[j_iter,1:mstar]+gam1*ones(1,mstar))')); endif; else; if m==1; F_CAST1=(YY[J_ITER,1]-X_MAT[J_ITER,1:p]*PHI)*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(st_mat[.,p+1:mstar+p]*rev((gam1*ones(1,mstar))')); i_iter=1; do until i_iter>p; F_CAST1=F_CAST1 + MU_MAT[.,p+1-i_iter]*PHI[i_iter] + lambda*(st_mat[.,p+1-i_iter:mstar+p-i_iter]*rev((gam1*ones(1,mstar))'))*PHI[i_iter]; i_iter=i_iter+1; endo; elseif m==2; F_CAST1=(YY[J_ITER,1]-X_MAT[J_ITER,1:p]*PHI)*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(1-st_mat[.,mstar+p+1]).*(st_mat[.,p+1:mstar+p]*rev((gam1*ones(1,mstar))')); i_iter=1; do until i_iter>p; F_CAST1=F_CAST1 + MU_MAT[.,p+1-i_iter]*PHI[i_iter] + lambda*(1-st_mat[.,mstar+1-i_iter]).*(st_mat[.,p+1-i_iter:mstar+p-i_iter]*rev((gam1*ones(1,mstar))'))*PHI[i_iter]; i_iter=i_iter+1; endo; elseif m==3; F_CAST1=(YY[J_ITER,1]-X_MAT[J_ITER,1:p]*PHI)*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(st_mat[.,p+1:mstar+p]*rev((x_mat[j_iter,1:mstar]+gam1*ones(1,mstar))')); i_iter=1; do until i_iter>p; F_CAST1=F_CAST1 + MU_MAT[.,p+1-i_iter]*PHI[i_iter] + lambda*(st_mat[.,p+1-i_iter:mstar+p-i_iter]*rev((x_mat[j_iter,1+i_iter:mstar+i_iter]+gam1*ones(1,mstar))'))*PHI[i_iter]; i_iter=i_iter+1; endo; endif; endif; VAR_L=VAR*ONES(DMNSION,1); PROB_DD=PR_TRF .* PROB__; @-----------------------------------------------------@ TMP=PROB_DD; i_iter=1; do until i_iter>mstar+p; tmp=tmp[1:ceil(0.5*rows(tmp))]+tmp[ceil(0.5*rows(tmp))+1:rows(tmp)]; i_iter=i_iter+1; endo; PR_STL0[J_ITER,1]=TMP[1,1]; @Pr[S_t=0|Y_t]@ @------------------------------------------------------@ if error_dist==1; PR_VL=(1./SQRT(2.*PI.*VAR)).*EXP(-0.5*F_CAST1.*F_CAST1./VAR).*PROB_DD; else; PR_VL=studt_d(f_cast1,0,VAR,NU_HAT).*PROB_DD; endif; @PR[S_t,S_{T-1},S_{T-2},...,S_{T-5},Y_t|Y_{t-1}]@ PR_VAL=SUMC(PR_VL); @f(y_t| Y_{t-1})@ LIK=-1*LN(PR_VAL); PRO_=PR_VL/PR_VAL; @Pr[S_t, S_{t-1},S_{t-2},...,S_{t-5} | Y_t]@ @-------------------------------------------------------@ TMP=PRO_; i_iter=1; do until i_iter>mstar+p; tmp=tmp[1:ceil(0.5*rows(tmp))]+tmp[ceil(0.5*rows(tmp))+1:rows(tmp)]; i_iter=i_iter+1; endo; PR_STT0[J_ITER,1]=TMP[1,1]; @Pr[S_t=0|Y_t]@ @-------------------------------------------------------@ PROB__T=PRO_[1:DMNSION/2,1]+PRO_[DMNSION/2+1:DMNSION,1]; @Pr[S_t, S_{t-1},....,S_{t-4}| Y_t]@ PROB__=VECR(PROB__T~PROB__T); J_ITER = J_ITER+1; endo; retp(PR_STT0, pr_stl0); endp; @++++++++++++++++++++++++++++++++++++++++++++++++++++++++@ proc get_probs(PRMTR1,q_stop); local prmtr, ppr,qpr,prob__0,prob__1,QQ, lik, j_iter, pr__0_l,pr__1_l, F_cast, var_L, pr_val,likv,phi,PSIX, vecp,st,st1,st2,st3,st4,ST5,ST6,ST7,ST8,ST9,ST10,ST11,ST12,ST13, pr_tr,pr_trf1,pr_trf,prob__t,prob__,pro_,pr_vl,j,psi1,psi2,var_c, delta0,DELTA1,MU0,MU1,st_k,st_k1,st_k2,st_k3,st_k4, f_cast1,f_cast2,PR_VL1,PR_VL2,pr_trf0, PR_TRF2,PR_TRF3,PR_TRF4,PR_TRF5,PR_TRF6,PR_TRF7,PR_TRF8,psic,psiL,PR_STL0,PR_STT0, TMPRY1,TMPRY2,SM_PRL,TMP_P0,SM_PR0,JJJ,MU_MAT,D_MAT,FLT_PRN, F1,F2,TMP00,TMP0,SM_PR00,prob_dd,VAR,TMP,OUT_MAT,A,EN,omga1,omga2, omga3,omga4,cdr1,cdr2,cdr3,cdr4,cdr1_m,cdr2_m,cdr3_m,cdr4_m,var1,var2, mu01,mu11,mu02,mu12,mu_add,pi_bar,mu_bar,d_mu,mu_pb,d_mu_pb,p00,p11,d_gam,lambda,gam0,gam1, mu_bar_pb,var_pb,d_gam_pb,i_iter,nu_hat,mu_bar_pb2; PRMTR=TRANS(PRMTR1); p00=prmtr[1]; @Pr[St=0/St-1=0]@ p11=prmtr[2]; @Pr[St=1/St-1=1]@ pi_bar=(1-p00)/(2-p11-p00); @Pr[St=1]@ var=prmtr[3]^2; mu_bar=prmtr[4]; d_gam=prmtr[5]; lambda=prmtr[6]; gam0=mu_bar; gam1=d_gam; @ recession, vs. boom@ if p==0; phi=0; else; phi= prmtr[7:6+p]; endif; @Structural break parameters@ d_gam_pb=d_gam; mu_bar_pb=mu_bar; mu_bar_pb2=mu_bar; var_pb=var; if sb==8; mu_bar_pb=prmtr[6+p+1]*mu_bar; mu_bar_pb2=mu_bar_pb; d_gam_pb=prmtr[6+p+2]*d_gam; var_pb=prmtr[6+p+2]^2*var; endif; if error_dist==2; if PRMTR[rows(prmtr)]<1000; NU_HAT=PRMTR[rows(prmtr)]; else; nu_hat=1000; endif; endif; MU_MAT = ST_MAT*gam1 + (ONES(DMNSION,NO_ST))*gam0; PR_TR=(p00~ (1-p11))| ((1-p00)~ p11); A = (eye(2)-pr_tr)|ones(1,2); EN=(0|0|1); if det(A'A) le 0; PROB__T=0.5|0.5; else; PROB__T = INV(A'A)*A'EN; @PR[S_t=0]|PR[S_t=1], 2x1 steady-state PROBABILITIES@ endif; PR_TRF0=VEC(PR_TR); i_iter=1; do until i_iter>mstar+p-1; PROB__T=VECR(PROB__T~PROB__T).*PR_TRF0; PR_TRF0=PR_TRF0|PR_TRF0; i_iter=i_iter+1; endo; PR_TRF =PR_TRF0; PROB__=VECR(PROB__T~PROB__T); PR_STT0=ZEROS(T,1); @WILL save Pr[S_t=0|Y_{t}@ PR_STL0=ZEROS(T,1); @WILL save Pr[S_t=0|Y_{t-1}@ LIKV=0.0; J_ITER=1; do until J_ITER>q_stop; if j_iter>breakdate1; gam0=mu_bar_pb; MU_MAT=ST_MAT*gam1 + (ONES(DMNSION,NO_ST))*gam0; endif; if j_iter>breakdate2; VAR=var_pb; gam0=mu_bar_pb2; gam1=d_gam_pb; MU_MAT=ST_MAT*gam1 + (ONES(DMNSION,NO_ST))*gam0; endif; if p==0; if m==1; F_CAST1=(YY[J_ITER,1])*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(st_mat[.,p+1:mstar+p]*rev((gam1*ones(1,mstar))')); elseif m==2; F_CAST1=(YY[J_ITER,1])*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(1-st_mat[.,mstar+p+1]).*(st_mat[.,p+1:mstar+p]*rev((gam1*ones(1,mstar))')); elseif m==3; F_CAST1=(YY[J_ITER,1])*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(st_mat[.,p+1:mstar+p]*rev((x_mat[j_iter,1:mstar]+gam1*ones(1,mstar))')); endif; else; if m==1; F_CAST1=(YY[J_ITER,1]-X_MAT[J_ITER,1:p]*PHI)*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(st_mat[.,p+1:mstar+p]*rev((gam1*ones(1,mstar))')); i_iter=1; do until i_iter>p; F_CAST1=F_CAST1 + MU_MAT[.,p+1-i_iter]*PHI[i_iter] + lambda*(st_mat[.,p+1-i_iter:mstar+p-i_iter]*rev((gam1*ones(1,mstar))'))*PHI[i_iter]; i_iter=i_iter+1; endo; elseif m==2; F_CAST1=(YY[J_ITER,1]-X_MAT[J_ITER,1:p]*PHI)*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(1-st_mat[.,mstar+p+1]).*(st_mat[.,p+1:mstar+p]*rev((gam1*ones(1,mstar))')); i_iter=1; do until i_iter>p; F_CAST1=F_CAST1 + MU_MAT[.,p+1-i_iter]*PHI[i_iter] + lambda*(1-st_mat[.,mstar+1-i_iter]).*(st_mat[.,p+1-i_iter:mstar+p-i_iter]*rev((gam1*ones(1,mstar))'))*PHI[i_iter]; i_iter=i_iter+1; endo; elseif m==3; F_CAST1=(YY[J_ITER,1]-X_MAT[J_ITER,1:p]*PHI)*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(st_mat[.,p+1:mstar+p]*rev((x_mat[j_iter,1:mstar]+gam1*ones(1,mstar))')); i_iter=1; do until i_iter>p; F_CAST1=F_CAST1 + MU_MAT[.,p+1-i_iter]*PHI[i_iter] + lambda*(st_mat[.,p+1-i_iter:mstar+p-i_iter]*rev((x_mat[j_iter,1+i_iter:mstar+i_iter]+gam1*ones(1,mstar))'))*PHI[i_iter]; i_iter=i_iter+1; endo; endif; endif; /* 2^6x1 */ VAR_L=VAR*ONES(DMNSION,1); PROB_DD=PR_TRF .* PROB__; @-----------------------------------------------------@ TMP=PROB_DD; i_iter=1; do until i_iter>mstar+p; tmp=tmp[1:ceil(0.5*rows(tmp))]+tmp[ceil(0.5*rows(tmp))+1:rows(tmp)]; i_iter=i_iter+1; endo; PR_STL0[J_ITER,1]=TMP[1,1]; @Pr[S_t=0|Y_t]@ @------------------------------------------------------@ if error_dist==1; PR_VL=(1./SQRT(2.*PI.*VAR)).*EXP(-0.5*F_CAST1.*F_CAST1./VAR).*PROB_DD; else; PR_VL=studt_d(f_cast1,0,VAR,NU_HAT).*PROB_DD; endif; @PR[S_t,S_{T-1},S_{T-2},...,S_{T-5},Y_t|Y_{t-1}]@ PR_VAL=SUMC(PR_VL); @f(y_t| Y_{t-1})@ LIK=-1*LN(PR_VAL); PRO_=PR_VL/PR_VAL; @Pr[S_t, S_{t-1},S_{t-2},...,S_{t-5} | Y_t]@ @-------------------------------------------------------@ TMP=PRO_; i_iter=1; do until i_iter>mstar+p; tmp=tmp[1:ceil(0.5*rows(tmp))]+tmp[ceil(0.5*rows(tmp))+1:rows(tmp)]; i_iter=i_iter+1; endo; PR_STT0[J_ITER,1]=TMP[1,1]; @Pr[S_t=0|Y_t]@ @-------------------------------------------------------@ PROB__T=PRO_[1:DMNSION/2,1]+PRO_[DMNSION/2+1:DMNSION,1]; @Pr[S_t, S_{t-1},....,S_{t-4}| Y_t]@ PROB__=VECR(PROB__T~PROB__T); J_ITER = J_ITER+1; endo; retp(pro_); endp; @======================================================================@ @================================================================@ @This procedure calculates the Pr[St|It], where It is info at time t@ proc smooth(prmtr1); local j,st1,st,t0,ystar, sig_v,st2,st3,st4,st5,st6,st7, st8,pr_tr,pr_trf1, pr_trf2,a,en,var_lm,PR_STMAT,TAU0,dtt0,stt0, var_l,pr_trf0,pr_trf, pr_trf3, pr_trf4, pr_trf5, pr_trf6, pr_trf7, f1,prob__1,prob__0,prob__t,prob__, jjj,prob_dd,pr_vl,pr_val,pro_,s_t,p0,p1,TMP,PROB__T0, phi_a,phi_b,intcpt, dum_star, gamma_vec, tstar, pt0_vec, pt0, pt1_vec, ptl_vec, ptt_vec, pro_all, pt1, f_cast0, f_cast1, pl_00, pl_01, pl_10, pl_11, pr_vl0, pr_vl1, prob_, data_tmp, mu_mat, pr_all9, pr_all10,PR_STP0,PR_STP1,PR_STP2,PR_STP3,PR_STP4,PR_STP5,PR_STP6,PR_STP7,PR_STPF, pt00,pt01,pt10,pt11,prmtr,ppr,qpr,var,mu0,mu1,omga1,phi,pr_stt0,pr_stt1,j_iter,pr_stl0,pr_sm0,pr_sm1, mu01,mu11,mu02,mu12,var1,var2,mu_add,pi_bar,mu_bar,d_mu,mu_pb,d_mu_pb,p00,p11,d_gam,lambda,gam0,gam1, mu_bar_pb,var_pb,d_gam_pb,i_iter,pr_all,nu_hat,mu_bar_pb2; pt0_vec = {}; pt1_vec = {}; ptl_vec = {}; ptt_vec = {}; PRMTR=TRANS(PRMTR1); p00=prmtr[1]; @Pr[St=0/St-1=0]@ p11=prmtr[2]; @Pr[St=1/St-1=1]@ pi_bar=(1-p00)/(2-p11-p00); @Pr[St=1]@ var=prmtr[3]^2; mu_bar=prmtr[4]; d_gam=prmtr[5]; lambda=prmtr[6]; gam0=mu_bar; gam1=d_gam; @ recession, vs. boom@ if p==0; phi=0; else; phi= prmtr[7:6+p]; endif; @Structural break parameters@ d_gam_pb=d_gam; mu_bar_pb=mu_bar; mu_bar_pb2=mu_bar; var_pb=var; if sb==8; mu_bar_pb=prmtr[6+p+1]*mu_bar; mu_bar_pb2=mu_bar_pb; d_gam_pb=prmtr[6+p+2]*d_gam; var_pb=prmtr[6+p+2]^2*var; endif; if error_dist==2; if PRMTR[rows(prmtr)]<1000; NU_HAT=PRMTR[rows(prmtr)]; else; nu_hat=1000; endif; endif; MU_MAT = ST_MAT*gam1 + (ONES(DMNSION,NO_ST))*gam0; PR_TR=(p00~ (1-p11))| ((1-p00)~ p11); A = (eye(2)-pr_tr)|ones(1,2); EN=(0|0|1); if det(A'A) le 0; PROB__T=0.5|0.5; else; PROB__T = INV(A'A)*A'EN; @PR[S_t=0]|PR[S_t=1], 2x1 steady-state PROBABILITIES@ endif; PR_TRF0=VEC(PR_TR); i_iter=1; do until i_iter>mstar+p-1; PROB__T=VECR(PROB__T~PROB__T).*PR_TRF0; PR_TRF0=PR_TRF0|PR_TRF0; i_iter=i_iter+1; endo; PR_TRF =PR_TRF0; PROB__=VECR(PROB__T~PROB__T); J_ITER=1; do until J_ITER>T; if j_iter>breakdate1; gam0=mu_bar_pb; MU_MAT=ST_MAT*gam1 + (ONES(DMNSION,NO_ST))*gam0; endif; if j_iter>breakdate2; VAR=var_pb; gam0=mu_bar_pb2; gam1=d_gam_pb; MU_MAT=ST_MAT*gam1 + (ONES(DMNSION,NO_ST))*gam0; endif; if p==0; if m==1; F_CAST1=(YY[J_ITER,1])*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(st_mat[.,p+1:mstar+p]*rev((gam1*ones(1,mstar))')); elseif m==2; F_CAST1=(YY[J_ITER,1])*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(1-st_mat[.,mstar+p+1]).*(st_mat[.,p+1:mstar+p]*rev((gam1*ones(1,mstar))')); elseif m==3; F_CAST1=(YY[J_ITER,1])*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(st_mat[.,p+1:mstar+p]*rev((x_mat[j_iter,1:mstar]+gam1*ones(1,mstar))')); endif; else; if m==1; F_CAST1=(YY[J_ITER,1]-X_MAT[J_ITER,1:p]*PHI)*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(st_mat[.,p+1:mstar+p]*rev((gam1*ones(1,mstar))')); i_iter=1; do until i_iter>p; F_CAST1=F_CAST1 + MU_MAT[.,p+1-i_iter]*PHI[i_iter] + lambda*(st_mat[.,p+1-i_iter:mstar+p-i_iter]*rev((gam1*ones(1,mstar))'))*PHI[i_iter]; i_iter=i_iter+1; endo; elseif m==2; F_CAST1=(YY[J_ITER,1]-X_MAT[J_ITER,1:p]*PHI)*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(1-st_mat[.,mstar+p+1]).*(st_mat[.,p+1:mstar+p]*rev((gam1*ones(1,mstar))')); i_iter=1; do until i_iter>p; F_CAST1=F_CAST1 + MU_MAT[.,p+1-i_iter]*PHI[i_iter] + lambda*(1-st_mat[.,mstar+1-i_iter]).*(st_mat[.,p+1-i_iter:mstar+p-i_iter]*rev((gam1*ones(1,mstar))'))*PHI[i_iter]; i_iter=i_iter+1; endo; elseif m==3; F_CAST1=(YY[J_ITER,1]-X_MAT[J_ITER,1:p]*PHI)*ONES(DMNSION,1)-MU_MAT[.,mstar+p+1]-lambda*(st_mat[.,p+1:mstar+p]*rev((x_mat[j_iter,1:mstar]+gam1*ones(1,mstar))')); i_iter=1; do until i_iter>p; F_CAST1=F_CAST1 + MU_MAT[.,p+1-i_iter]*PHI[i_iter] + lambda*(st_mat[.,p+1-i_iter:mstar+p-i_iter]*rev((x_mat[j_iter,1+i_iter:mstar+i_iter]+gam1*ones(1,mstar))'))*PHI[i_iter]; i_iter=i_iter+1; endo; endif; endif; VAR_L=(VAR)*ONES(DMNSION,1); PROB_DD=PR_TRF .* PROB__; @Pr[St-8,...,St-1,St | I_t-1]@ @DMNSIONx1@ if error_dist==1; PR_VL=(1./SQRT(2.*PI.*VAR)).*EXP(-0.5*F_CAST1.*F_CAST1./VAR).*PROB_DD; else; PR_VL=studt_d(f_cast1,0,VAR,NU_HAT).*PROB_DD; endif; PR_VAL=SUMC(PR_VL); PRO_ = PR_VL/PR_VAL; @Pr[St-8,...,St-1,St | I_t]@ @DMNSIONx1@ PROB__T = PRO_[1:DMNSION/2,1] + PRO_[DMNSION/2+1:DMNSION,1]; @Pr[St-7,...,St-1,St | I_t]@ PROB__=VECR(PROB__T~PROB__T); TMP=PRO_; i_iter=1; do until i_iter>mstar+p; tmp=tmp[1:ceil(0.5*rows(tmp))]+tmp[ceil(0.5*rows(tmp))+1:rows(tmp)]; i_iter=i_iter+1; endo; pt0 = tmp[1]; @Pr[St=0|I_t]@ @1x1@ pt1 = tmp[2]; @Pr[St=1|I_t]@ @1x1@ ptl_vec=ptl_vec|prob_dd'; @Pr[St-8,...,St-1,St | I_t-1]@ @TxDMNSION@ ptt_vec=ptt_vec|pro_'; @Pr[St-8,...,St-1,St | I_t]@ @TxDMNSION@ j_iter = j_iter+1; endo; pr_all=pro_; @At time T, set Pr[St-7,...,St,St+1 | I_T] for next iteration when t=T-1@ @DMNSIONx1@ PR_STP0=0|1|0|1; i_iter=1; do until i_iter>mstar+p; PR_STP0=PR_STP0|PR_STP0; i_iter=i_iter+1; endo; PR_STPF=PR_STP0; @Pr[S_T|I_T]@ pr_sm0=pt0; pr_sm1=pt1; jjj = T - 1; do until jjj<1; @Smoother@ pr_all = (pr_all|pr_all).*vecr(ptt_vec[jjj,.]'~ptt_vec[jjj,.]').*(PR_TRF|PR_TRF)./(ptl_vec[jjj+1,.]'|ptl_vec[jjj+1,.]'); @Pr[St-8,...,St,St+1|I_T]=Pr[St-7,...,St,St+1|I_T]*(Pr[St-8,...,St|I_t]*Pr[St+1|St])/Pr[St-7,...,St,St+1|I_t for jjj=T-1@ @Pr[St-8,...,St,St+1|I_T,St+2]=Pr[St-7,...,St,St+1|I_T,St+2]*(Pr[St-8,...,St|I_t]*Pr[St+1|St])/Pr[St-7,...,St,St+1|I_t for jjjmstar+p; tmp=tmp[1:ceil(0.5*rows(tmp))]+tmp[ceil(0.5*rows(tmp))+1:rows(tmp)]; i_iter=i_iter+1; endo; @Pr[St|I_T]@ pt0 = tmp[1]; pt1 = tmp[2]; pr_sm0=pt0|pr_sm0; pr_sm1=pt1|pr_sm1; jjj=jjj-1; endo; retp(pr_sm0); @This proc returns Pr[S_t=0|Y_T]@ endp; /*********************************************************/ /*********************************************************/ proc TRANS(c0); @ constraining values of reg. coeff.@ local c1; c1=c0; c1[1:2,.]=exp(c0[1:2,.])./ (1+exp(c0[1:2,.])); @p00 and p11@ c1[3]=exp(c0[3]); @var@ c1[5]=-0.1-exp(c0[5]); @d_gamma@ @c1[6]=-exp(c0[6]);@ @lambda@ if sb==8; c1[6+p+1]=exp(c0[6+p+1]); @d_gamma_pb@ c1[6+p+2]=exp(c0[6+p+2]); @var_pb@ endif; if error_dist==2; c1[rows(c0)]=1+10000*(exp(c0[rows(c0)])/(1+exp(c0[rows(c0)]))); endif; retp(c1); endp; @========================================================================@