/* =============================================================== File: a:\SRTVP.opt MODEL: Time Varying Parameter model for stock return predictability. Measurement Equation: r(t,t+1) - MU = beta_k(t) * (r(t-k,t) - k * MU) + e_t e_t ~ N(0, RR) State (Transition) Equation: beta_k(t) = beta_k(t-1) + v_t v_t ~ N(0, QQ) ===> 3 parameters (MU, RR, and QQ) to be estimated WRITTEN BY JAMES MORLEY, DEPT. OF ECONOMICS, UNIVERSITY OF WASHINGTON MODIFIED FROM A PROGRAM BY CHANG_JIN KIM, DEPT. OF ECONOMICS, KOREA UNIVERSITY ================================================================ */ @ Version 5/9/98@ @ Used for PHASE 6@ New; library optmum,pgraph; #include optmum.ext; #include pgraph.ext; optset; graphset; bfixed=0; k=96; @Select holding period (in months) for k-period returns@ r=2; @Select equal (r=1) or value (r=2) weighted returns@ @Select decile(j-2) (r=j)@ load data[852,12]=c:\CRSPD96.DAT; /* column 1: CRSP monthly index of NYSE equal-weighted total real return 2: CRSP monthly index of NYSE value-weighted total real return 3-12: decile returns Sample period -- 1926:1--1986:12,... 732 observations */ T=rows(data); @Create the Data@ @----------------------------------------------------------------@ @Equal Weighted one month and lagged k-period returns@ ew=ln(1+data[1:T,1]); ew_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; ew_kl[j+1,1]=sumc(ew[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Value Weighted one month and lagged k-period returns@ vw=ln(1+data[1:T,2]); vw_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; vw_kl[j+1,1]=sumc(vw[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 1 one month and lagged k-period returns@ d1=ln(1+data[1:T,3]); d1_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d1_kl[j+1,1]=sumc(d1[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 2 one month and lagged k-period returns@ d2=ln(1+data[1:T,4]); d2_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d2_kl[j+1,1]=sumc(d2[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 3 one month and lagged k-period returns@ d3=ln(1+data[1:T,5]); d3_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d3_kl[j+1,1]=sumc(d3[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 4 one month and lagged k-period returns@ d4=ln(1+data[1:T,6]); d4_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d4_kl[j+1,1]=sumc(d4[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 5 one month and lagged k-period returns@ d5=ln(1+data[1:T,7]); d5_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d5_kl[j+1,1]=sumc(d5[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 6 one month and lagged k-period returns@ d6=ln(1+data[1:T,8]); d6_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d6_kl[j+1,1]=sumc(d6[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 7 one month and lagged k-period returns@ d7=ln(1+data[1:T,9]); d7_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d7_kl[j+1,1]=sumc(d7[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 8 one month and lagged k-period returns@ d8=ln(1+data[1:T,10]); d8_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d8_kl[j+1,1]=sumc(d8[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 9 one month and lagged k-period returns@ d9=ln(1+data[1:T,11]); d9_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d9_kl[j+1,1]=sumc(d9[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @Decile 10 one month and lagged k-period returns@ d10=ln(1+data[1:T,12]); d10_kl=zeros(T,1); @create space for lagged k-period returns@ j=k; do until j>T-1; d10_kl[j+1,1]=sumc(d10[j-k+1:j,1]); @k-period returns are lagged here@ j=j+1; endo; @------------------------------------------------------------------@ format /m1 /rd 8,6; output file=c:\p6\r2k96.out reset; output off; data=ew~vw~d1~d2~d3~d4~d5~d6~d7~d8~d9~d10~ew_kl~vw_kl~ d1_kl~d2_kl~d3_kl~d4_kl~d5_kl~d6_kl~d7_kl~d8_kl~d9_kl~d10_kl; data=data[k+1:T,1:24]; S=T-k; @Subsample after adjusting for lagged X@ if r<2; yy=data[.,1]; xx=data[.,13]; @ew returns@ goto skip1; endif; if r<3; yy=data[.,2]; xx=data[.,14]; @vw returns@ goto skip1; endif; if r<4; yy=data[.,3]; xx=data[.,15]; @vw returns@ goto skip1; endif; if r<5; yy=data[.,4]; xx=data[.,16]; @vw returns@ goto skip1; endif; if r<6; yy=data[.,5]; xx=data[.,17]; @vw returns@ goto skip1; endif; if r<7; yy=data[.,6]; xx=data[.,18]; @vw returns@ goto skip1; endif; if r<8; yy=data[.,7]; xx=data[.,19]; @vw returns@ goto skip1; endif; if r<9; yy=data[.,8]; xx=data[.,20]; @vw returns@ goto skip1; endif; if r<10; yy=data[.,9]; xx=data[.,21]; @vw returns@ goto skip1; endif; if r<11; yy=data[.,10]; xx=data[.,22]; @vw returns@ goto skip1; endif; if r<12; yy=data[.,11]; xx=data[.,23]; @vw returns@ goto skip1; endif; if r<13; yy=data[.,12]; xx=data[.,24]; @vw returns@ goto skip1; endif; skip1: @================= Initialize Global Variables============@ START=25; @Ignore first START-1 observations in calculating lik @ PRMTR_NO=3; @ >>Number of parameters to be estimated. >>@ TVP_NO=1; @ >>Number of Time-varying Coefficients >>@ pr_1=3.098578 ;pr_2=0.4;pr_3=0.006413 ; @for r=2,k=120, scale is different, pr_2=0.7@ @initial values of parameters@ PRMTR_IN=pr_1|pr_2|pr_3; @OUTPUT@ output on; if r<2; "Equal Weighted Real Returns"; goto skip2; endif; if r<3; "Value Weighted Real Returns"; goto skip2; endif; if r<4; "Decile(1) Real Returns"; goto skip2; endif; if r<5; "Decile(2) Real Returns"; goto skip2; endif; if r<6; "Decile(3) Real Returns"; goto skip2; endif; if r<7; "Decile(4) Real Returns"; goto skip2; endif; if r<8; "Decile(5) Real Returns"; goto skip2; endif; if r<9; "Decile(6) Real Returns"; goto skip2; endif; if r<10; "Decile(7) Real Returns"; goto skip2; endif; if r<11; "Decile(8) Real Returns"; goto skip2; endif; if r<12; "Decile(9) Real Returns"; goto skip2; endif; if r<13; "Decile(10) Real Returns"; goto skip2; endif; skip2: " "; "==MODEL======================================================="; " Measurement (Observation) Equation:"; " r(t,t+1) - MU = beta_k(t) * (r(t-k,t) - k * MU) + e_t"; " e_t ~ N(0, RR)"; " "; " State (Transition) Equation:"; " beta_k(t) = beta_k(t-1) + v_t"; " v_t ~ N(0, QQ)"; "==============================================================="; " "; "==Estimation Details==========================================="; "k=";;k; "Start=";;START;;" => t=1 @ January, ";;1926+(k+Start-1)/12; " (first 12 observations are used to start up Kalman filter)"; "==============================================================="; " "; output off; @ Maximum Likelihood Estimation @ {xout,fout,gout,cout}=optmum(&lik_fcn,PRMTR_IN); @ prmtr estimates, -log lik value, Grandient, Hessian@ "Calculating Hessian..... Please be patient!!!!"; hout0=hessp(&lik_fcn,xout); hout=inv(hout0); PRM_FNL=TRANS(xout); @ Estimated coefficients, constrained@ 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; "==FINAL OUTPUT================================================="; "likelihood value is ";; -fout; "Estimated parameters are:"; "Sigma_e Sigma_v Mu"; prm_fnl';"(";;xout';;")"; "Standard errors of parameters are:"; sd_fnl'; "==============================================================="; " "; output off; {F_SS,tvp_coef}=SAVE_OUT(xout); OUTPUT FILE=c:\p6\r2k96.DTA RESET; F_SS~tvp_coef; OUTPUT OFF; begwind; window(2,2,0); xlabel("month"); xindex=seqa(1,1,rows(f_ss)); setwind(1); title("forecast error"); xy(xindex, f_ss[.,1]); nextwind; title("conditional variance"); xy(xindex, f_ss[.,2]); nextwind; title("beta_t|t"); xy(xindex,tvp_coef[.,1]~tvp_coef[.,1]+2*tvp_coef[.,2]^0.5~ zeros(rows(f_ss),1)~tvp_coef[.,1]-2*tvp_coef[.,2]^0.5~zeros(rows(f_ss),1)); nextwind; title("beta_t|T"); xy(xindex,tvp_coef[.,3]~tvp_coef[.,3]+2*tvp_coef[.,4]^0.5~ zeros(rows(f_ss),1)~tvp_coef[.,3]-2*tvp_coef[.,4]^0.5~zeros(rows(f_ss),1)); endwind; begwind; xlabel("month"); title("Comparison of Beta's (fixed, filtered, and smoothed)"); xy(xindex, tvp_coef[.,1]~ tvp_coef[.,3]~ tvp_coef[.,1]~ bfixed*ones(rows(f_ss),1)~ zeros(rows(f_ss),1)~ zeros(rows(f_ss),1)~ zeros(rows(f_ss),1)); endwind; @Standardized Forecast Errors@ sfe=f_ss[.,3]./f_ss[.,4]^0.5; @(v_t|T)/(f_t|T)@ @Quasi-Standardized Forecast Errors@ qsfe=f_ss[.,3]./prm_fnl[1]; @(v_t|T)/(sig_e)@ begwind; _pbartyp = { 1 7}; title("Histogram of Standardized Forecast Errors"); xlabel("v_t|T/f_t|T"); call hist(sfe,25); endwind; @Descriptive Statistics@ mean_sfe=meanc(sfe); med_sfe=median(sfe); max_sfe=maxc(sfe); min_sfe=minc(sfe); std_sfe=stdc(sfe); skew_sfe=(1/rows(sfe))*sumc(sfe.^3)/((1/rows(sfe))*sumc(sfe.^2))^3/2; kurto_sfe=(1/rows(sfe))*sumc(sfe.^4)/((1/rows(sfe))*sumc(sfe.^2))^2; w_sfe=rows(sfe)*(((skew_sfe^2)/6) + (((kurto_sfe-3)^2)/24)); pval_sfe = cdfchic(w_sfe,2); mean_qsfe=meanc(qsfe); med_qsfe=median(qsfe); max_qsfe=maxc(qsfe); min_qsfe=minc(qsfe); std_qsfe=stdc(qsfe); skew_qsfe=(1/rows(qsfe))*sumc(qsfe.^3)/((1/rows(qsfe))*sumc(qsfe.^2))^3/2; kurto_qsfe=(1/rows(qsfe))*sumc(qsfe.^4)/((1/rows(qsfe))*sumc(qsfe.^2))^2; w_qsfe=rows(qsfe)*(((skew_qsfe^2)/6) + (((kurto_qsfe-3)^2)/24)); pval_qsfe = cdfchic(w_qsfe,2); output file=c:\p6\r2k96.out; output on; "==Descriptive Statistics======================================="; "Standardized Forecast Errors, (v_t|T)/(f_t|T): "; " Mean = ";;mean_sfe; " Median = ";;med_sfe; " Maximum = ";;max_sfe; " Minimum = ";;min_sfe; " Standard Deviation = ";;std_sfe; " Skewness = ";;skew_sfe; " Kurtosis = ";;kurto_sfe; " "; " Jarque-Bera = ";;w_sfe; " p-value = ";;pval_sfe; " "; "Quasi-Standardized Forecast Errors. (v_t|T)/(sig_e): "; " Mean = ";;mean_qsfe; " Median = ";;med_qsfe; " Maximum = ";;max_qsfe; " Minimum = ";;min_qsfe; " Standard Deviation = ";;std_qsfe; " Skewness = ";;skew_qsfe; " Kurtosis = ";;kurto_qsfe; " "; " Jarque-Bera = ";;w_qsfe; " p-value = ";;pval_qsfe; "==============================================================="; output off; @END OF MAIN PROGRAM@ @========================================================================@ @========================================================================@ PROC LIK_FCN(PRMTR1); local BETA_TT, P_TT, QQ, LIK, J_ITER, XT, BETA_TL, P_TL, SS,K_GAIN, F_CAST, F, BETA_LL, P_LL,RR,YT,LIK_MAT,PRMTR,MU; PRMTR=TRANS(PRMTR1); Locate 15,1; PRMTR'; RR=PRMTR[1,1]^2; QQ=PRMTR[2,1]^2; MU=PRMTR[3,1]; F=1; BETA_LL=0; @initial guess for coefficients@ P_LL=EYE(TVP_NO)*10; @uncertainty associated with initial guess@ LIK_MAT=ZEROS(S,1); J_ITER = 1; DO UNTIL J_ITER>S; XT=XX[J_ITER,.] - k*MU; YT=YY[J_ITER,.] - MU; BETA_TL = F * BETA_LL; P_TL = F * P_LL * F' +QQ; F_CAST = YT - XT * BETA_TL; SS = XT * P_TL * XT' + RR; BETA_TT = BETA_TL + (P_TL*XT'/SS) * F_CAST; P_TT = (EYE(TVP_NO) - (P_TL*XT'/SS)* XT ) * P_TL; BETA_LL=BETA_TT; P_LL=P_TT; IF J_ITER < START; goto skip; endif; LIK_MAT[J_ITER,1] = 0.5*(LN(2*pi*SS) + (F_CAST^2)/SS); skip: J_ITER = J_ITER+1; ENDO; RETP(SUMC(LIK_MAT[START:S,1])); ENDP; @========================================================================@ PROC (2) = SAVE_OUT(PRMTR1); local BETA_TT, P_TT, QQ, LIK, J_ITER, XT, BETA_TL, P_TL, SS,K_GAIN, F_CAST, F,DTA_MAT,DTA_MAT2,BETA_LL, P_LL,SIG_E,SIG_1,PRMTR, RR,YT,MU,BTT_MAT,PTL_MAT,PTT_MAT,BETA_PALL,P_PALL,BETA_TALL, P_TALL,SS_ALL,F_CASTALL; PRMTR=TRANS(PRMTR1); LOCATE 15,1; PRMTR'; BTT_MAT=zeros(S,1); PTL_MAT=zeros(S,1); PTT_MAT=zeros(S,1); DTA_MAT=zeros(S,4); DTA_MAT2=zeros(S,4); RR=PRMTR[1,1]^2; QQ=PRMTR[2,1]^2; MU=PRMTR[3,1]; F=1; BETA_LL=0; @initial guess for coefficients@ P_LL=EYE(TVP_NO)*10; @uncertainty associated with initial guess@ J_ITER = 1; DO UNTIL J_ITER>S; XT=XX[J_ITER,.] - k*MU; YT=YY[J_ITER,.] - MU; BETA_TL = F * BETA_LL; P_TL = F * P_LL * F' +QQ; F_CAST = YT - XT * BETA_TL; SS = XT * P_TL * XT' + RR; BETA_TT = BETA_TL + (P_TL*XT'/SS) * F_CAST; P_TT = (EYE(TVP_NO) - (P_TL*XT'/SS) * XT ) * P_TL; BTT_MAT[J_ITER,.]=BETA_TT; PTL_MAT[J_ITER,.]=P_TL; PTT_MAT[J_ITER,.]=P_TT; DTA_MAT[J_ITER,1:2]=F_CAST~SS; @Forecast error and conditional variance@ DTA_MAT2[J_ITER,1:2]=BETA_TT~P_TT; BETA_LL=BETA_TT; P_LL=P_TT; J_ITER = J_ITER+1; ENDO; BETA_PALL=BETA_TT; P_PALL=P_TT; DTA_MAT2[S,3:4]=BETA_TT~P_TT; J_ITER = S-1; DO UNTIL J_ITER