/* Morley (2007) JMCB UC model of Income and Consumption with cointegration and correlated innvoations Replicates Table 3 */ new; library optmum, pgraph; format /m1 /rd 9,5; load data_im[214,4]=yc.txt; @ gdp_pc dpy_pc cnds_pc cnds_pc_apx 1952:Q1 to 2005:Q @ data=data_im[9:214,1]~data_im[9:214,3]; y=100*ln(data); T=rows(y); @ Maximum Likelihood Estimation @ START=2; prior=100; PRMTR_IN=zeros(11,1)|(y[1,1]-y[1,2])/100|y[1,2]/100|1; prmtr_in={ 0.49144 -1.19353 0.76793 -6.54179 0.874 1.51205 1.97532 1.50067 0.00914 -0.01329 -0.00769 0.5 1.00079 }; prmtr_in={ 1 -4 0.74382 -5.07080 0.51159 -0.25900 0.41104 7.31927 1.02369 -1.50885 -0.76931 -0.16347 0.96781 }; prmtr_in=prmtr_in'; output file=uc_yc.out reset; output on; "Starting values:"; prmtr_in'; output off; @Initial values of the parameters @ {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_fnl =sqrt(diag(cov)); @Standard errors of the estimated coefficients@ sd_out=sqrt(diag(cov0)); output on; "==FINAL OUTPUT================================================"; "likelihood value is ";; -fout; "code";;cout; " "; "Estimated parameters are:"; ""; prm_fnl'|sd_fnl'; "=============================================================="; "Pre-transformed estimates are:"; "(";;XOUT';;")"; OUTPUT OFF; {DATA,FORCST} = filter(xout); output file=uc_yc.dta reset; output on; data[.,1]~data[.,3]~data[.,5]~forcst[.,1:2]; output off; phi1_1=prm_fnl[1]; phi2_1=prm_fnl[2]; phi1_2=prm_fnl[3]; phi2_2=prm_fnl[4]; f_1=phi1_1~phi2_1|1~0; f_2=phi1_2~phi2_2|1~0; irf_fnl={}; irf=1; psi_ll=0; psi_l=1; j=1; do until j>40; psi_t=phi1_1*psi_l+phi2_1*psi_ll; irf=irf|psi_t; psi_ll=psi_l; psi_l=psi_t; j=j+1; endo; irf_fnl=irf_fnl~irf; irf=1; psi_ll=0; psi_l=1; j=1; do until j>40; psi_t=phi1_2*psi_l+phi2_2*psi_ll; irf=irf|psi_t; psi_ll=psi_l; psi_l=psi_t; j=j+1; endo; irf_fnl=irf_fnl~irf; hlp=0.5*ones(rows(irf_fnl),1); @half lives@ hlm=-0.5*ones(rows(irf_fnl),1); output file=model1_eig.out reset; output on; eig(f_1); abs(eig(f_1)); eig(f_2); abs(eig(f_2)); output off; output file=ucy_yc_irf.dta reset; output on; irf_fnl; output off; begwind; window(3,2,0); xlabel("month"); setwind(1); xlabel("month"); title("income"); xy(SEQA(1,1,ROWS(Y[start:T,1])),Y[start:T,1]); nextwind; xlabel("month"); title("consumption"); xy(SEQA(1,1,ROWS(Y[start:T,1])),Y[start:T,2]); nextwind; xlabel("month"); title("Permanent Income"); _pdate = ""; XY(SEQA(1,1,ROWS(DATA)),DATA[.,5]); nextwind; xlabel("month"); title("Permanent Consumption"); _pdate = ""; XY(SEQA(1,1,ROWS(DATA)),DATA[.,5]+prm_fnl[12]); nextwind xlabel("month"); title("Transitory Income"); _pdate = ""; XY(SEQA(1,1,ROWS(DATA)),DATA[.,1]~zeros(rows(data),1)); nextwind xlabel("month"); title("Transitory Consumption"); _pdate = ""; XY(SEQA(1,1,ROWS(DATA)),DATA[.,3]~zeros(rows(data),1)); endwind; begwind; window(2,1,0); xlabel("quarter"); setwind(1); xlabel("quarter"); title("irf"); xy(SEQA(1,1,ROWS(irf_fnl)),irf_fnl[.,1]~hlp~zeros(rows(irf_fnl),1)~hlm); nextwind; xlabel("quarter"); title("irf"); xy(SEQA(1,1,ROWS(irf_fnl)),irf_fnl[.,2]~hlp~zeros(rows(irf_fnl),1)~hlm); endwind; end; @ END OF MAIN PROGRAM @ @========================================================================@ @========================================================================@ PROC lik_fcn(PRMTR1); LOCAL F,FSTAR,H,BETA_LL,P_LL,BETA_TL,P_TL,BETA_TT,P_TT, VECP_LL,vt,ft, VAL,LIK_MAT,J_ITER, prmtr, phi11,phi12,phi21,phi22,ddd,sig11,sig22,sigvv,sig12,sig1v, sig2v,y0,c0,Q,QSTAR,mu,a,beta_c,aaa; PRMTR=TRANS(PRMTR1); @LOCATE 20,1; PRMTR';@ phi11=prmtr[1]; phi12=prmtr[2]; phi21=prmtr[3]; phi22=prmtr[4]; ddd=prmtr[5]; sig11=prmtr[6]^2; sig22=prmtr[7]^2; sigvv=prmtr[8]^2; sig12=prmtr[9]*sqrt(sig11*sig22); sig1v=prmtr[10]*sqrt(sig11*sigvv); sig2v=prmtr[11]*sqrt(sig22*sigvv); aaa=prmtr[12]; beta_c=prmtr[13]; F= PHI11~PHI12~0~0~0| 1~0~0~0~0| 0~0~PHI21~PHI22~0| 0~0~1~0~0| 0~0~0~0~1; FSTAR= PHI11~PHI12~0~0| 1~0~0~0| 0~0~PHI21~PHI22| 0~0~1~0; MU=(0~0~0~0~ddd)'; H= 1~0~0~0~1| 0~0~1~0~beta_c; Q= sig11~0~sig12~0~sig1v| 0~0~0~0~0| sig12~0~sig22~0~sig2v| 0~0~0~0~0| sig1v~0~sig2v~0~sigvv; QSTAR= sig11~0~sig12~0| 0~0~0~0| sig12~0~sig22~0| 0~0~0~0; A=(0~aaa)'; BETA_LL=(0~0~0~0~947.5)'; VECP_LL=INV(EYE(16) - FSTAR.*.FSTAR)*vec(QSTAR); P_LL= VECP_LL[1,1]~0~VECP_LL[3,1]~0~0| 0~0~0~0~0| VECP_LL[9,1]~0~VECP_LL[11,1]~0~0| 0~0~0~0~0| 0~0~0~0~prior; LIK_MAT=ZEROS(T,1); J_ITER = 1; DO UNTIL J_ITER>T; BETA_TL = MU + F * BETA_LL ; P_TL = F * P_LL * F' + Q ; vt=y[j_iter,1:2]' - H * BETA_TL - A; @prediction error@ ft= H * P_TL * H' ; @variance of forecast error@ @invft=inv(ft);@ 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)^2)*det(ft)) +0.5*vt'*inv(ft)*vt; BETA_LL=BETA_TT; P_LL=P_TT; J_ITER = J_ITER+1; ENDO; VAL = SUMC(LIK_MAT[START:T]); RETP(VAL); ENDP; @>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@ PROC (2)=filter(PRMTR1); LOCAL F,FSTAR,H,BETA_LL,P_LL,BETA_TL,P_TL,BETA_TT,P_TT, VECP_LL,vt,ft, VAL,LIK_MAT,J_ITER, prmtr, DDD1,DDD2,DDD3,DDD4,DDD5,invft,FE_MAT,BETA_MAT,MU,A,a12,a21,a31,a41, PHI11,PHI12,PHI13,FCST_MAT,a11,a22,a33,a44,a55, PHI21,PHI22,PHI23,a42, G,GSTAR,a51,a52,a43,a53,a54, PHI31,PHI32,PHI41,PHI42,PHI51,PHI52,Q,QSTAR, phi11,phi12,phi21,phi22,ddd,sig11,sig22,sigvv,sig12,sig1v, sig2v,y0,c0,mu,a,beta_c,aaa; PRMTR=TRANS(PRMTR1); LOCATE 20,1; PRMTR'; phi11=prmtr[1]; phi12=prmtr[2]; phi21=prmtr[3]; phi22=prmtr[4]; ddd=prmtr[5]; sig11=prmtr[6]^2; sig22=prmtr[7]^2; sigvv=prmtr[8]^2; sig12=prmtr[9]*sqrt(sig11*sig22); sig1v=prmtr[10]*sqrt(sig11*sigvv); sig2v=prmtr[11]*sqrt(sig22*sigvv); aaa=prmtr[12]; beta_c=prmtr[13]; F= PHI11~PHI12~0~0~0| 1~0~0~0~0| 0~0~PHI21~PHI22~0| 0~0~1~0~0| 0~0~0~0~1; FSTAR= PHI11~PHI12~0~0| 1~0~0~0| 0~0~PHI21~PHI22| 0~0~1~0; MU=(0~0~0~0~ddd)'; H= 1~0~0~0~1| 0~0~1~0~beta_c; Q= sig11~0~sig12~0~sig1v| 0~0~0~0~0| sig12~0~sig22~0~sig2v| 0~0~0~0~0| sig1v~0~sig2v~0~sigvv; QSTAR= sig11~0~sig12~0| 0~0~0~0| sig12~0~sig22~0| 0~0~0~0; A=(0~aaa)'; BETA_LL=(0~0~0~0~947.5)'; VECP_LL=INV(EYE(16) - FSTAR.*.FSTAR)*vec(QSTAR); P_LL= VECP_LL[1,1]~0~VECP_LL[3,1]~0~0| 0~0~0~0~0| VECP_LL[9,1]~0~VECP_LL[11,1]~0~0| 0~0~0~0~0| 0~0~0~0~prior; BETA_MAT=ZEROS(T,5); FCST_MAT=ZEROS(T,2); J_ITER = 1; DO UNTIL J_ITER>T; BETA_TL = MU + F * BETA_LL ; P_TL = F * P_LL * F' + Q; vt=y[j_iter,1:2]' - H * BETA_TL - A ; @prediction error@ ft= H * P_TL * H'; @variance of forecast error@ @invft=inv(ft);@ BETA_TT= BETA_TL + P_TL * H' * inv(ft) * vt; P_TT= P_TL - P_TL * H' * inv(ft) * H * P_TL; BETA_LL=BETA_TT; P_LL=P_TT; FCST_MAT[J_ITER,.]=vt'; BETA_MAT[J_ITER,.]=BETA_TT'; J_ITER = J_ITER+1; ENDO; RETP(BETA_MAT[START:T,.],FCST_MAT[START:T,.]); ENDP; @================================================================@ PROC trans(c0); @ constraining values of reg. coeff.@ local c1,aaa,ccc,ggg,qqq,c11,c22,c33,c44,c55,c31,c42,c51, c21,c31,c32; c1=c0; c1[12]=100*c0[12]; c11=exp(-c0[6]); c22=exp(-c0[7]); c33=exp(-c0[8]); c21=c0[9]; c31=c0[10]; c32=c0[11]; ggg= c11~0~0| c21~c22~0| c31~c32~c33; qqq=ggg*ggg'; c1[6]=sqrt(qqq[1,1]); c1[7]=sqrt(qqq[2,2]); c1[8]=sqrt(qqq[3,3]); c1[9]=qqq[2,1]/sqrt(qqq[1,1]*qqq[2,2]); c1[10]=qqq[3,1]/sqrt(qqq[1,1]*qqq[3,3]); c1[11]=qqq[3,2]/sqrt(qqq[2,2]*qqq[3,3]); aaa=c0[1]./(1+abs(c0[1])); ccc=(1-abs(aaa))*c0[2]./(1+abs(c0[2]))+abs(aaa)-aaa^2; c1[1]=2*aaa; c1[2]= -1* (aaa^2+ccc) ; aaa=c0[3]./(1+abs(c0[3])); ccc=(1-abs(aaa))*c0[4]./(1+abs(c0[4]))+abs(aaa)-aaa^2; c1[3]=2*aaa; c1[4]= -1* (aaa^2+ccc) ; retp(c1); endp;