/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ uc_ur.opt This program generates Table 4 in Morley, Nelson, and Zivot, "Why Are Beveridge-Nelson and Unobserved-Component Decompositions of GDP So Different?" (forthcoming, RESTAT). If you use this program, please cite the paper. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ new; library optmum, pgraph; format /m1 /rd 9,6; @Real GDP 1947.1 - 1998.2; From CITIBASE@ load data[206,1]=lngdpq.txt; @Real GDP 1947.1 - 1998.2; From CIT@ yy= 100*data; dyy=yy[2:195]-yy[1:194]; t=rows(yy); output file=uc_ur.out reset; output off; @ Maximum Likelihood Estimation @ START=2; @ startup values for VEVD of likelihood @ @Initial values of the parameters @ PRMTR_IN={ 0 0 0 0 0 0 }; @PRMTR_IN={ 1.2368 0.7487 -0.8392 0.8155 1.3418 -0.7058 };@ @ use above starting values with no constraints to eliminate imaginary standard errors@ @ Note: numerical standard errors for GAUSS differ from computer to computer@ PRMTR_IN=PRMTR_IN'; @ maximize likelihood function @ {xout,fout,gout,cout}=optmum(&lik_fcn,PRMTR_IN); @ prmtr estimates, -log lik value, Grandient, code@ prm_fnl=trans(xout); "Calculating Hessian..... Please be patient!!!!"; hessn0=hessp(&lik_fcn,xout); cov0=inv(hessn0); grdn_fnl=gradfd(&TRANS,xout); cov=grdn_fnl*cov0*grdn_fnl'; SD =sqrt(diag(cov)); @Standard errors of the estimated coefficients@ output on; "==FINAL OUTPUT========================================================"; "initial values of prmtr is"; trans(prmtr_in)'; "=============================================================="; "likelihood value is ";; -fout; "code";;cout; "Estimated parameters are:"; prm_fnl'; XOUT'; "Standard errors of parameters are:"; sd'; OUTPUT OFF; f=prm_fnl[5]~prm_fnl[6]|1~0; ef=eig(f); output on; ef; output off; DATA = filter(xout); output file=uc_ur.dta reset; DATA; output off; XY(SEQA(1,1,ROWS(DATA)),YY[START:T,1]~DATA[.,1]); XY(SEQA(1,1,ROWS(DATA)),DATA[.,2]); @ END OF MAIN PROGRAM @ @========================================================================@ @========================================================================@ PROC LIK_FCN(PRMTR1); LOCAL SN,SE,SNE,F,F1,H,Q,Q1,R,BETA_LL,P_LL,BETA_TL,P_TL,BETA_TT,P_TT, P_LL1, vt,ft, VAL,LIK_MAT,J_ITER, phi1,phi2,vecpll,prmtr,muvec,mu; @ transform hyperparameters to impose contraints @ PRMTR=TRANS(PRMTR1); @LOCATE 10,1; PRMTR';@ SN=PRMTR[1]; @s.e. of the random walk component@ SE=PRMTR[2]; @s.e. of the AR component@ SNE=PRMTR[3]; @ COV(n,e) @ mu=PRMTR[4]; phi1=PRMTR[5]; phi2=PRMTR[6]; muvec = mu|0|0; @ drift vector @ F=(1~0~0)| @ transition matrix @ (0~phi1~phi2)| (0~1~0); F1 = submat(F,2~3,2~3); @ pick out I(0) part @ H=1~1~0; @ measurement equation @ Q= ZEROS(3,3); @ E[V(t)V(t)'] @ Q[1,1]=SN^2; Q[2,2]=SE^2; Q[1,2]=SNE; Q[2,1]=SNE; Q1 = submat(Q,2~3,2~3); @ cov matrix of I(0) part @ R= 0; @ cov matrix of error in measurement eqn @ BETA_LL=yy[1]|0|0; @ starting values @ @ variance matrix of initial state vector @ P_LL = zeros(3,3); P_LL1 = inv( (eye(4)-F1.*.F1) )*vec(Q1); P_LL1 = reshape(P_LL1',2,2); P_LL[1,1]=100000000; P_LL[2:3,2:3]=P_LL1; LIK_MAT=ZEROS(T,1); J_ITER = 1; DO UNTIL J_ITER>T; BETA_TL = muvec + F*BETA_LL ; P_TL = F*P_LL*F' + Q; vt=yy[j_iter,1] - H * BETA_TL ; @prediction error@ ft= H * P_TL * H' + R; @variance of forecast error@ BETA_TT= BETA_TL + P_TL * H' * inv(ft) * vt; P_TT= P_TL - P_TL * H' * inv(ft) * H * P_TL; LIK_MAT[J_ITER,1] = 0.5*(LN(2*pi*ft) + vt^2/ft); BETA_LL=BETA_TT; P_LL=P_TT; J_ITER = J_ITER+1; ENDO; VAL = SUMC(LIK_MAT[START:T]); RETP(VAL); ENDP; @>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@ PROC FILTER(PRMTR1); LOCAL SN,SE,SNE,F,F1,H,Q,Q1,R,BETA_LL,P_LL,BETA_TL,P_TL,BETA_TT,P_TT, P_LL1, vt,ft, VAL,LIK_MAT,J_ITER, phi1,phi2,vecpll,prmtr,muvec,mu, BETA_MAT; BETA_MAT=ZEROS(T,2); LIK_MAT=ZEROS(T,1); PRMTR=TRANS(PRMTR1); LOCATE 15,1; PRMTR'; SN=PRMTR[1]; @s.e. of the random walk component@ SE=PRMTR[2]; @s.e. of the AR component@ SNE=PRMTR[3]; @ COV(n,e) @ mu=PRMTR[4]; phi1=PRMTR[5]; phi2=PRMTR[6]; muvec = mu|0|0; @ drift vector @ F=(1~0~0)| @ transition matrix @ (0~phi1~phi2)| (0~1~0); F1 = submat(F,2~3,2~3); @ pick out I(0) part @ H=1~1~0; @ measurement equation @ Q= ZEROS(3,3); @ E[V(t)V(t)'] @ Q[1,1]=SN^2; Q[2,2]=SE^2; Q[1,2]=SNE; Q[2,1]=SNE; Q1 = submat(Q,2~3,2~3); @ cov matrix of I(0) part @ R= 0; @ cov matrix of error in measurement eqn @ BETA_LL=yy[1]|0|0; @ starting values @ @ variance matrix of initial state vector @ P_LL = zeros(3,3); P_LL1 = inv( (eye(4)-F1.*.F1) )*vec(Q1); P_LL1 = reshape(P_LL1',2,2); P_LL[1,1]=100000000; P_LL[2:3,2:3]=P_LL1; J_ITER = 1; DO UNTIL J_ITER>T; BETA_TL = muvec + F * BETA_LL ; P_TL = F * P_LL * F' + Q; vt=yy[j_iter,1] - H * BETA_TL ; @prediction error@ ft= H * P_TL * H' + R; @variance of forecast error@ BETA_TT= BETA_TL + P_TL * H' * inv(ft) * vt; P_TT= P_TL - P_TL * H' * inv(ft) * H * P_TL; BETA_MAT[J_ITER,.]=BETA_TT[1,1]~BETA_TT[2,1]; BETA_LL=BETA_TT; P_LL=P_TT; J_ITER = J_ITER+1; ENDO; RETP(BETA_MAT[START:T,.]); ENDP; @================================================================@ PROC TRANS(c0); @ constraining values of reg. coeff.@ /* c0 = (sn,se,sne,mu,ph1,ph2) constraints imposed 1. cov(n_t,e_t) is pd 2. AR(2) is stationary */ local c1,p,varcov,aaa,ccc,e1,e2; c1 = c0; c1[1:2]=exp(-c0[1:2]); /* constrain cov(n_t,e_t) to be pd using Cholesky factorization */ p = zeros(2,2); p[1,1] = exp(-c0[1]); p[2,1] = c0[3]; p[2,2] = exp(-c0[2]); varcov = p'p; c1[1]=sqrt(varcov[1,1]); c1[2]=sqrt(varcov[2,2]); c1[3]=varcov[2,1]; aaa=c0[5]./(1+abs(c0[5])); ccc=(1-abs(aaa))*c0[6]./(1+abs(c0[6]))+abs(aaa)-aaa^2; c1[5]=2*aaa; c1[6]= -1* (aaa^2+ccc) ; retp(c1); endp;