/*********************************************************************************** This is the program file for Figures 1 and 2 in "Identifying Business Cycle Turning Points in Real Time" by Marcelle Chauvet and Jeremy Piger in the March/April 2003 issue of the Federal Reserve Bank of St. Louis Review. This is the GAUSS program files to reproduce figures 1 and 2 in the paper. I used windows GAUSS version 3.6 to run this program. This program also needs the optmum.lib module. Figures 1A and 2A use quarterly GDP data. The default setting for this program is to compute figure 2A. Before you can run the program, you need to download the real time quartly data set(rtqtrly.zip) from Philadelphia FED's web site: http://www.phil.frb.org/econ/forecast/readow.html Then you extract all the .wk1 files into a folder, and set "dir" (variable in this program) to the path of that folder. This program will automatically loop over all the spreadsheets and extract GDP data series from the files. If you want to reproduce figure 1A, keep data_type at 2, and set multi equal to 0. If you want to reproduce figures 1B and 2B, you need to modify the monthly employment data file (0303mcd.txt). See the header in emp.txt for instructions to set up the employment data files. You also need to change line 162, 165, and 168 to include the path to the directory where you have the employment data files. To recreate figure 1B, set data_type=1 and multi=0. To recreate figure 2B, set data_type=1 and multi=1. ***********************************************************************************/ /*======================================================================== MODEL: HAMILTON'S(1989,ECONOMETRICA) MARKOV-SWITCHING MODEL ***** THIS PROGRAM CONTAINS KIM'S(1994) SMOOTHING ALGORITHM***** (y_t-mu_{s_t})= phi_1*(y_{t-1}-mu_{s_{t-1}})+... +phi_4*(y_{t-4}-mu_{s_{t-4}}) + e_t e_t -- i.i.d. N(0,sigma^2) mu_{st} = (1-S_t)*mu_0 + S_t* mu_1 Pr[S_t=1|S_{t-1}=1]=p Pr[S_t=0|S_{t-1}=0]=q WRITTEN BY CHANG-JIN KIM, based on Hamilton (1989) DEPT. OF ECONOMICS KOREA UNIVERSITY SEOUL, 136-701, KOREA cjkim@korea.ac.kr Modified by John Zhu and Jeremy Piger Research Department St. Louis Federal Reserve Bank =========================================================================*/ new; cls; library optmum,PGRAPH; /******************************** Before executing the program, please set "dir" to the directory with all the extracted spreedsheet files.(or where all the data files are stored) ********************************/ dir="H:\\piger\\review\\data\\"; data_type=2; @1=employment, 2=gdp@ multi=1; @multi=0: loop is off@ /**************************************** This part of the GAUSS program file reads in all the spreadsheet files in the zip file, and it extracts the second(GDP) column from each file in rtqrtly.zip and put them together into a single EXCEL file called gdp.xls. gdp.xls is an input file for the GAUSS program for the paper, "Identifying Business Cycle Turning Points in Real Time" by Marcelle Chauvet and Jeremy Piger in the March/April 2003 issue of the Federal Reserve Bank of St. Louis Review. Before executing the program, please set "dir" to the directory with all the extracted spreedsheet files. *****************************************/ if data_type .eq 2; "loading GDP spreadsheets. Please be patient..."; ftype=".wk1"; month="feb"|"may"|"aug"|"nov"; vnames={}; /*loops over all the wk1 files*/ year=1965; endyear=2002; month_ind=4; /*starts in nov*/ max_rows=400; column=3; /*gdp is on column 3*/ data=zeros(max_rows,1); data=miss(data,0); do while year <= endyear; if year%100 < 10; fmat="0%lf"; else; fmat="%lf"; endif; yr=ftos(year%100,fmat,1,0); do while month_ind%5 .ne 0; fname=dir$+month[month_ind%5]$+yr$+ftype; screen off; {x,names}=import(fname,0,0); screen on; data=data~data[.,cols(data)]; x=x[.,column]; i=1; do while i<=rows(x); if (x[i]<-1000000); x[i]=-1000000; endif; i=i+1; endo; data[1:rows(x),cols(data)-1]=x; vnames=vnames|month[month_ind%5]$+yr; month_ind=month_ind+1; endo; month_ind=month_ind+1; year=year+1; endo; data=miss(data,-1000000); data=data[.,1:cols(data)-1]; screen off; export(data,dir$+"gdp.xls",vnames); screen on; "loading complete."; endif; /******************************************************************* This program reads in three files from the data folder: emp64_74.txt, emp75_89.txt, and emp90.txt. version switching: This program can either run on the aug 2002 data serie or the entire real time data set. set multi=1 if you want to loop for the entire real time data. set multi=0 if you just want the single aug 2002 serie. The real time gdp data is avaliable at Philadelphia FED's website: http://www.phil.frb.org/econ/forecast/readow.html ********************************************************************/ if multi .eq 1; rt={}; /*this will save the last element of pr_tt0 after every iteration*/ endif; if data_type .eq 1; @1=employment@ "loading montly employment data..."; data={}; load data[] = H:\\piger\\review\\data\\emp64_74.txt; @monthly, U.S.employment@ col=data[1]; data=data[2:rows(data)]; @number of columns in the file@ emp=reshape(data,rows(data)/col,col); load data[]= H:\\piger\\review\\data\\emp75_89.txt; @monthly, U.S.employment@ col=data[1]; data=data[2:rows(data)]; @number of columns in the file@ emp2=reshape(data,rows(data)/col,col); load data[]= H:\\piger\\review\\data\\emp90.txt; @monthly, U.S.employment@ col=data[1]; data=data[2:rows(data)]; @number of columns in the file@ emp3=reshape(data,rows(data)/col,col); @delete the date column@ emp=emp[.,2:cols(emp)]; emp2=emp2[.,2:cols(emp2)]; emp3=emp3[.,2:cols(emp3)]; @merge the three data matrices@ data=zeros(rows(emp3),cols(emp)); data=miss(data,0); data[1:rows(emp),1:cols(emp)]=emp; emp=data; data=zeros(rows(emp3),cols(emp2)); data=miss(data,0); data[1:rows(emp2),1:cols(emp2)]=emp2; emp2=data; data=emp~emp2~emp3; "loading complete."; elseif data_type .eq 2; @2=gdp@ "reloading GDP data..."; screen off; {data,names}=import(dir$+"gdp.xls",0,0); emp=data; screen on; "GDP data loading complete."; else; "invalid data_type."; end; endif; if multi .eq 0; @test a single serie: aug 2002@ if data_type .eq 1; data=data[.,cols(data)-3]; elseif data_type .eq 2; data=data[.,cols(data)-1]; else; "invalid data type."; end; endif; data=packr(data); @delete the missing values@ else; emp=data; /*emp now contains the real time data series*/ endif; @================= Initialize Global Variables======1959.1-89.9======@ START=1; @1952:4.....@ mu0_s = 1; mu1_s = -.8; sig1_s = 1; p00_s = 3; p11_s = 2; PRMTR_IN=p11_s~p00_s~sig1_s~mu0_s~mu1_s; PRMTR_IN=PRMTR_IN'; i=1; flag=1; do while (i <= cols(emp)) .and (flag .eq 1); if multi .eq 1; data=packr(emp[.,i]); else; flag=0; /*if it is single serie, exit after one iteration*/ endif; @for employment data, we start at 1947.11, and gdp data at 1947.q4@ if data_type .eq 1; @monthly employment data@ yy_d=(ln(data[98:rows(data)])-ln(data[97:rows(data)-1]))*100; elseif data_type .eq 2; @quarterly gdp data@ yy_d=(ln(data[3:rows(data)])-ln(data[2:rows(data)-1]))*100; else; "invalid data type"; end; endif; t0=rows(yy_d); LAG_AR=4; NO_ST=lag_ar+1; @ NUMBER OF STATES TO BE CONSIDERED@ DMNSION=2^NO_ST; st_mat=zeros(DMNSION,NO_ST); j=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,.]=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; yy=yy_d[lag_ar+1:t0,1]; x_mat= yy_d[Lag_ar:T0-1,1] ~ yy_d[lag_ar-1:T0-2,1]~ yy_d[lag_ar-2:T0-3,1]~ yy_d[lag_ar-3:T0-4,1]; T=rows(yy); @ Maximum Likelihood Estimation @ @==================================================@ {xout,fout,gout,cout}=optmum(&lik_fcn,PRMTR_in); PRM_FNL=TRANS(xout); @Estimated coefficients, constrained@ "==FINAL OUTPUT========================================================"; "initial values of prmtr is"; trans(prmtr_in)'; "=============================================================="; "code is--------";;cout; "likelihood value is ";; -1*fout; "Estimated parameters are:"; prm_fnl'; " "; xout'; output off; "Calculating Hessian..... Please be patient!!!!"; hout0=hessp(&lik_fcn,xout); hout=inv(hout0); grdn_fnl=gradfd(&TRANS,xout); Hsn_fnl=grdn_fnl*hout*grdn_fnl'; SD_fnl =sqrt(diag(Hsn_fnl)); @Standard errors of the estimated coefficients@ output on; "Standard errors of parameters are:"; sd_fnl'; "==============================================================="; output off; {pr_tt0,pr_tl0}=FILTER(XOUT); @Pr[S_t=0|Y_t] and Pr[S_t=0|Y_{t-1}]@ smooth0=SMOOTH(pr_tt0,pr_tl0); @Pr[S_t=0|Y_T], Smoothed probabilities@ FLTR=(1-pr_tt0)~(1-pr_tl0)~(1-smooth0); if multi .eq 0; xy(seqa(1,1,rows(pr_tt0)),(1-pr_tt0)); endif; @xy(seqa(1,1,rows(smooth0)),(1-smooth0));@ if multi .eq 1; rt=rt|(1-pr_tt0[rows(pr_tt0)]); endif; if i > 2; xy(seqa(1,1,rows(rt)),(rt)); endif; prmtr_in=xout; i=i+1; endo; end; @ END OF MAIN PROGRAM @ @========================================================================@ @========================================================================@ PROC LIK_FCN(PRMTR1); local prmtr, ppr,qpr,prob__0,prob__1,QQ, lik, j_iter, pr__0_l,pr__1_l, F_cast, var_L,pr_vl, pr_val,likv,phi,PSIC,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_trf7,pr_trf0, PR_TRF2,PR_TRF3,PR_TRF4,PR_TRF5,PR_TRF6,psic,psiL, TMPRY1,TMPRY2,SM_PRL,TMP_P0,SM_PR0,JJJ,MU_MAT,D_MAT,FLT_PRN, F1,F2,TMP00,TMP0,SM_PR0,SM_PR00,prob_dd,VAR,A,EN; PRMTR=TRANS(PRMTR1); LOCATE 15,1; PRMTR'; /* DEFINE PARAMETERS */ PPR=PRMTR[1,1]; @Pr[St=1/St-1=1]@ QPR=PRMTR[2,1]; @Pr[St=0/St-1=0]@ PHI=0|0|0|0; VAR=PRMTR[3,1]^2; MU0=PRMTR[4,1]; @ recession, vs. boom@ MU1=PRMTR[5,1]; @ recession, vs. boom@ MU_MAT=ST_MAT*MU1 + (ONES(DMNSION,NO_ST)-ST_MAT)*MU0; /* A MATRIX OF TRANSITION PROBABILITIES */ PR_TR=(QPR~ (1-PPR))| ((1-QPR)~ PPR); /* INITIALIZING THE FILTER WITH STEADY-STATE PROBABILITIES */ A = (eye(2)-pr_tr)|ones(1,2); EN=(0|0|1); if rows(packr(a)) .ne rows(a); retp(999999); endif; PROB__T = INV(A'A)*A'EN; @PR[S_t=0]|PR[S_t=1], 2x1 steady-state PROBABILITIES@ PR_TRF0=VEC(PR_TR); PR_TRF1=PR_TRF0|PR_TRF0; PR_TRF2=PR_TRF1|PR_TRF1; PR_TRF =PR_TRF2|PR_TRF2; PROB__T=VECR(PROB__T~PROB__T).*PR_TRF0; /*PR[S_{-2},S_{-3}|I_0) 4x1*/ PROB__T=VECR(PROB__T~PROB__T).*PR_TRF1; /*PR[S_{-1},S_{-2},S_{-3}|I_0) 8x1*/ PROB__T=VECR(PROB__T~PROB__T).*PR_TRF2; /*PR[S_{0},S_{-1},S_{-2},S_{-3}|I_0) 16x1 */ PROB__=VECR(PROB__T~PROB__T); /*2^5 x 1*/ LIKV=0.0; J_ITER=1; DO UNTIL J_ITER>T; F_CAST1=(YY[J_ITER,1]-X_MAT[J_ITER,.]*PHI)*ONES(DMNSION,1) -(MU_MAT[.,5] - MU_MAT[.,4]*PHI[1,1] - MU_MAT[.,3]*PHI[2,1] - MU_MAT[.,2]*PHI[3,1] - MU_MAT[.,1]*PHI[4,1] ); /* 2^5x1 */ VAR_L=VAR*ONES(DMNSION,1); /* 2^5x1 */ PROB_DD=PR_TRF .* PROB__; /* Pr[S_t,S_{t-1},S_{t-2},S_{t-3},S_{t-4} | I(t-1)]*/ PR_VL=(1./SQRT(2.*PI.*VAR_L)).*EXP(-0.5*F_CAST1.*F_CAST1./VAR_L).*PROB_DD; /* 2^5x1 */ /* Joint density of y_t,S_t,..,S_{t-4} 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^5 conditional densities */ LIK=-1*LN(PR_VAL); PRO_=PR_VL/PR_VAL; /* Pr[S_t,S_{t-1},S_{t-2},S_{t-3},S_{t-4} | 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-2},S_{t-3}| Y_t] */ /* 2^4x1*/ PROB__=VECR(PROB__T~PROB__T); /* 2^5x1 */ LIKV = LIKV+LIK; J_ITER = J_ITER+1; ENDO; LOCATE 2,35;"LIKV=";;LIKV; 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_vl, pr_val,likv,phi,PSIC,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_trf7,pr_trf0, PR_TRF2,PR_TRF3,PR_TRF4,PR_TRF5,PR_TRF6,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_PR0,SM_PR00,prob_dd,VAR,TMP,OUT_MAT,A,EN; PRMTR=TRANS(PRMTR1); LOCATE 15,1; PRMTR'; PPR=PRMTR[1,1]; @Pr[St=1/St-1=1]@ QPR=PRMTR[2,1]; @Pr[St=0/St-1=0]@ PHI=0|0|0|0; VAR=PRMTR[3,1]^2; MU0=PRMTR[4,1]; @ recession, vs. boom@ MU1=PRMTR[5,1]; @ recession, vs. boom@ PR_TR=(QPR~ (1-PPR))| ((1-QPR)~ PPR); MU_MAT=ST_MAT*MU1 + (ONES(DMNSION,NO_ST)-ST_MAT)*MU0; /* FOR UNCONDITIONAL PROBABILITIES */ 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 UNCONDITIONAL PROBABILITIES@ PR_TRF0=VEC(PR_TR); PR_TRF1=PR_TRF0|PR_TRF0; PR_TRF2=PR_TRF1|PR_TRF1; PR_TRF =PR_TRF2|PR_TRF2; PROB__T=VECR(PROB__T~PROB__T).*PR_TRF0; @2^2 x 1@ PROB__T=VECR(PROB__T~PROB__T).*PR_TRF1; @2^3 x 1@ PROB__T=VECR(PROB__T~PROB__T).*PR_TRF2; @2^4 x 1@ PROB__=VECR(PROB__T~PROB__T); @2^5 x 1@ 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; F_CAST1=(YY[J_ITER,1]-X_MAT[J_ITER,.]*PHI)*ONES(DMNSION,1) -(MU_MAT[.,5] - MU_MAT[.,4]*PHI[1,1] - MU_MAT[.,3]*PHI[2,1] - MU_MAT[.,2]*PHI[3,1] - MU_MAT[.,1]*PHI[4,1] ); VAR_L=VAR*ONES(DMNSION,1); PROB_DD=PR_TRF .* PROB__; @-----------------------------------------------------@ TMP=PROB_DD; TMP=TMP[1:16]+TMP[17:32]; TMP=TMP[1:8]+TMP[9:16]; TMP=TMP[1:4]+TMP[5:8]; TMP=TMP[1:2]+TMP[3:4]; PR_STL0[J_ITER,1]=TMP[1,1]; @Pr[S_t=0|Y_t]@ @------------------------------------------------------@ PR_VL=(1./SQRT(2.*PI.*VAR_L)).*EXP(-0.5*F_CAST1.*F_CAST1./VAR_L).*PROB_DD; @PR[S_t,S_{T-1},S_{T-2},S_{T-3},S_{T-4},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-3},S_{t-4} | Y_t]@ @-------------------------------------------------------@ TMP=PRO_; TMP=TMP[1:16]+TMP[17:32]; TMP=TMP[1:8]+TMP[9:16]; TMP=TMP[1:4]+TMP[5:8]; TMP=TMP[1:2]+TMP[3:4]; 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-2},S_{t-3}| Y_t]@ PROB__=VECR(PROB__T~PROB__T); J_ITER = J_ITER+1; ENDO; RETP(PR_STT0,PR_STL0); ENDP; @======================================================================@ @======================================================================@ PROC SMOOTH(pr_tt0,pr_tl0); @pr_TT0 contains Pr[S_t|Y_t]@ @pr_TL0 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_tt1,pr_tl1; PPR=PRM_fnl[1,1]; @Pr[St=1/St-1=1]@ QPR=PRM_fnl[2,1]; @Pr[St=0/St-1=0]@ pr_tt1=1-pr_tt0; pr_tl1=1-pr_tl0; pr_sm0=pr_tt0; @ pr_sm0 will contain Pr[S_t|Y_T]@ pr_sm1=pr_tt1; j_iter=T-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_sm0); @This proc returns Pr[S_t=0|Y_T]@ endp; @++++++++++++++++++++++++++++++++++++++++++++++++++++++++@ PROC TRANS(c0); @ constraining values of reg. coeff.@ local c1,m,u,d1,d2,d3,d4,d5,d6; c1=c0; c1[1:2,.]=exp(c0[1:2,.])./ (1+exp(c0[1:2,.])); retp(c1); endp;