/* GIBBS SAMPLING APPROACH TO HAMILTON'S (1989) AR(0) MODEL WITH A MARKOV-SWITCHING MEAN, APPLIED TO REAL GNP. BASED ON THE ORIGINAL CODE WRITTEN BY ALTERT AND CHIB (JBES, 1993). ------------------------------------------------------------------- A MULTI-MOVE GIBBS SAMPLING IS INCORPORATED. */ new; #lineson; library pgraph; format /m1 /rd 9,5; load data[195,1]=gdp4795.prn; @Quarter, U.S. real GDP@ @1947.1 -- 1995.3 @ ytt=(ln(data[21:195,1])-ln(data[20:194,1]))*100; @21:152; 1952.1 -- 1984.4 [hamilton data set]@ t=rows(ytt); @================= PRIOR PARAMETERS ====================@ U1_01_=2; U1_00_=8; @FOR A1TT@ U1_10_=1; U1_11_=9; @FOR B1TT@ R0_M=EYE(2)*1; @For MU0, MU1@ T0_M=0|1; D0_=0; V0_=0; @FOR SIG_0^2@ @========INITIAL VALUES FOR THE GIBBS SAMPLER ===========@ A1TT=1-0.8; B1TT=1-0.8; MU0TT=0; MU1TT=0.5; SIG2TT=1; @===== DESIGN FOR THE GIBBS SAMPLER ===========@ N0 = 2000; @NUMBER OF DRAWS TO LEAVE OUT@ MM = 10000; @NUMBER OF DRAWS@ LL = 5; @EVERY L-TH VALUE IS USED@ CAPN = N0 + MM; @TOTAL NUMBER OF DRAWS@ @===== STORAGE SPACES===========================@ MU0MM={}; MU1MM={}; SIG2MM={}; PMM={};QMM={}; sttmm=zeros(t,1); ITR=1; DO WHILE ITR LE CAPN; ITR;;"th iteration......."; @=================START SAMPLING================@ STT=GEN_ST; {MU0TT,MU1TT}=GEN_MU; SIG2TT=GEN_SIG; tranmat = switchg(stt+1,1|2); A1TT = rndb(tranmat[1,2]+u1_01_,tranmat[1,1]+u1_00_); B1TT = rndb(tranmat[2,1]+u1_10_,tranmat[2,2]+u1_11_); PTT=1-B1TT; QTT=1-A1TT; @===============================================@ locate 10,1; print "this is tranmat" tranmat; "mu";; mu0tt;;mu1tt; "SIG2";; SIG2TT; "P AND Q";; PTT;;QTT; IF ITR GT N0; SIG2MM=SIG2MM~SIG2TT; MU0MM=MU0MM~MU0TT; MU1MM=MU1MM~MU1TT; PMM=PMM~PTT; QMM=QMM~QTT; STTMM=STTMM+STT; ENDIF; ITR=ITR+1; ENDO; STTMM=STTMM/MM; FNL_MAT=SIG2MM'~MU0MM'~MU1MM'~PMM'~QMM'; SAVE FNL_MAT,STTMM; output file=gibs_ms0.dta reset; sttmm; output off; @======== calculation of posterior expectations and std deviations@ 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=GIBS_MS0.OUT RESET; "FOR SIG2, MU0, MU1, P, Q"; "================================================"; SORT_OUT; "PR[S_t=1|I(T)]"; "================================================"; STTMM; OUTPUT OFF; END; @========================================================@ @ A PROCEDURE THAT GENERATES MU0, MU1 @ @========================================================@ PROC (2) = GEN_MU; local ystar,xstar,v,c,accept,coef,root,rootmod, tstar,YSTARR,PHI_G,MU_G,PHI,MU; YSTAR=YTT; XSTAR=ONES(T,1)~STT; V = invpd(INV(R0_m) + (SIG2TT)^(-1)*XSTAR'XSTAR); MU = V*(INV(R0_m)*t0_m + (SIG2TT)^(-1)*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(MU_G[1],MU_G[2]); ENDP; @=========================================================================@ @==== PROCEDURE THAT GENERATES SIG2 =====================================@ @=========================================================================@ PROC GEN_SIG; local ystarr,ystar,xstar,v,psi,c,accept,psi_g,root,rootmod,nn,d,t2, sig2_v_g,coef,tstar,SIG2,Y_TMP; Y_TMP=YTT-MU0TT*ONES(T,1)-MU1TT*STT; TSTAR=ROWS(Y_tmp); nn = tstar + v0_; d = d0_ + (Y_tmp)'(Y_tmp); c = rndc(nn); t2 = c/d; SIG2=1/t2; RETP(SIG2); 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,st3,st4,phi,pr_tr,pr_trf1, pr_trf2,a,en, 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,TMP,PROB__T0; YSTAR=YTT; TSTAR=ROWS(YSTAR); FLT_PR=ZEROS(2,TSTAR); DMNSION=4; @For common component@ PPR=1-b1tt; @Pr[St=1/St-1=1]@ QPR=1-a1tt; @Pr[St=0/St-1=0]@ mu_0=mu0tt; mu_1=mu1tt; sig_v=SIG2TT; PR_TR=(QPR~ (1-PPR))| ((1-QPR)~ PPR); VAR_L=SIG_V*ONES(4,1); mu_mat=mu_0|(mu_0+mu_1)|mu_0|(mu_0+mu_1); @<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@ PR_TRF=VEC(PR_TR); A = (eye(2)-pr_tr)|ones(1,2); EN=(0|0|1); PROB__T = INV(A'A)*A'EN; @PR[S_t=0]|PR[S_t=1], 2x1 steady-state PROBABILITIES@ PROB__T0=PROB__T; PROB__=VECR(PROB__T~PROB__T); @++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++@ @BEGIN FILTERING ++++++++++++++++++++++++++++++++++++++++++++++++@ @++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++@ J_ITER=1; DO UNTIL J_ITER>TSTAR; f_cast1=ystar[j_iter,1]*ones(DMNSION,1)- mu_mat; 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); FLT_PR[.,J_ITER]=PROB__T; @Pr[St|Yt], 2x1@ J_ITER = J_ITER+1; ENDO; @=============== GENERATE S_T =====================================@ 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(S_T); 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;