/***************************************************************************/ /* */ /* Title: OBSFACT.G */ /* includes OBSFACT (procedure) */ /* OLD version before spr 95 modifications. Saved in case of disaster */ /* */ /* Written: John H. Cochrane */ /* by: University of Chicago Graduate School of Business */ /* 1101 E. 58th Street, Rosenwald 2xx */ /* Chicago, Ill 60637 */ /* ph#: (312) 702-3059 */ /* e-mail: fac_jhc@gsbacd.uchicago.edu */ /* */ /* Purpose: Does observable factor models with out searching for the */ /* technology and/or consumption parameters. Used to test */ /* whether CAPM, CRR factors price returns */ /* */ /* Data: Quarterly averaged returns created by redodat.gpg and GAUSS */ /* formatted (indicated by .fmt extension) */ /* */ /* */ /***************************************************************************/ /* ------------------------------------------------------------------------*/ /* PROC: obsfact.g */ /* */ /* USAGE: */ /* {jt,dof,neg_m} = obsfact(rassetx,factrs,xtrfact,inst,nwlags,ergrf,prnt, */ /* win,factlbl,xtrlbl,cfact,momcode,errbars); */ /* INPUTS: */ /* rassetx: excess returns to test */ /* factrs: observable factors */ /* xtrfact: extra factors, test whether not priced. to not use, use 0 */ /* inst: instruments. 0 of none. */ /* Make 47:1 to 90:4, program cuts to shape. */ /* nwlags: vector of newey west lags to try e.g. 0|4 */ /* ergrf: 3x1 vector (1,1)=do graph or E(r) vs. pred E(r)? */ /* (2,1)=view 1st stage? */ /* (3,1)=view 2nd stage? */ /* prnt: 1 or 0 print out estimates? */ /* win: 0: do regular 2 stage gmm. */ /* 1: do iterated gmm */ /* matrix: do second stage only, using given w matrix */ /* factlbl: character matrix of labels for factors */ /* xtrlbl: character matrix of labels for xtra factors */ /* cfact: 1 include scaled factors; 0 only use unscaled factors */ /* momcode: 1 include scaled returns; 0 only unscaled returns */ /* errbars: 0 no std error bars; 1 insert std error bars */ /* */ /* OUTPUTS: */ /* jt: JT matrix (see paper) */ /* dof: degrees of freedom */ /* neg_m: computed % of negative m's */ /* */ /* */ /* NOTE: if you use inst make sure inst is 47:1 to 90:2 and returns, */ /* factors are 47:3 to 90:4. */ /* ------------------------------------------------------------------------*/ /* ------------------------------------- */ /* Part I: (1) Call GAUSS library files */ /* (2) Procedure definition */ /* ------------------------------------- */ library pgraph,optmum,gmm; proc(6) = obsfact(rassetx,factrs,xtrfact,inst,nwlags,ergrf,prnt,win,factlbl, xtrlbl,cfact,momcode,errbars,preder2,erx2,grphcode, w_matrix); /* warning: can't use clearg with extra factors because of self-call */ local nf,nx,Er,C,b,b2,m,f,lgT,jT,W1,varmat,se,b1,m1,f1,gt1,jt1,doff,pval, preder,erx,T,W,j,k,label,s0,jtvec,allr,table,mask,fmt,blank,chi2, allfact,allfact1,ni,select,dofvec,trsh,stderr,tneg_m,neg_m,junk, ejunk,pjunk,gjunk,gt1_,iter,oldchi,test1,test2,lastw; /* ------------------------------------------------- */ /* Part II: (1) set up factors (demeaning & scaling) */ /* (2) set up returns (scale if required) */ /* ------------------------------------------------- */ tneg_m = 0; if ergrf == 0; ergrf = 0|0; endif; factrs = factrs-meanc(factrs)'; /* de mean factors */ nf = cols(factrs); /* nf = n. factors, NOT scaled*/ nx = 0; /* nx = n. extra, NOT scaled */ ni = 0; /* n of instruments NOT inc const*/ if inst /= 0; ni = cols(inst); endif; T = rows(rassetx); allfact = factrs; if xtrfact /= 0 ; /* set up extra factors */ xtrfact = xtrfact-meanc(xtrfact)'; nx = cols(xtrfact); allfact = factrs~xtrfact; endif; /* if wanted, scale factors */ if (inst /= 0) and (cfact == 1); allfact1 = allfact; j = 1; do while j <= cols(inst); allfact = allfact ~ ((allfact1 .* inst[.,j])-meanc(allfact1 .* inst[.,j])'); j = j+1; endo; endif; allr = rassetx; /* scale rets */ if momcode; j = 1; do while j <= cols(inst); allr = allr ~ (rassetx.*inst[.,j]); j = j+1; endo; endif; if cols(allr) < cols(allfact); "Unidentified system: more factors than returns"; retp(0,0,0,0,0); endif; /* -------------------------- */ /* Part III: First stage GMM */ /* -------------------------- */ /* ---------------- */ /* 1. Set up JT, gT */ /* ---------------- */ Er = meanc(allr); C = 1/T*(allr-Er')'allfact; /* Cov(R,f') */ if (win == 0) or (win == 1); /* dont' bother with first stage if given w */ W = eye(cols(allr)); b = - inv(C'W*C)*C'W*Er ; m = 1 + allfact*b; neg_m = counts(m,0)/rows(m); tneg_m = tneg_m + neg_m; f = (m .* allr); lgT = meanc(f); JT = T*lgT'W*lgT; /* -------------------------------------------- */ /* 2. Print first stage estimates, if requested */ /* -------------------------------------------- */ if prnt; "1. First stage GMM (identity weighting matrix)"; "Parameter estimates b"; let blank = " "; label = factlbl'; if xtrfact /= 0; label = factlbl'~xtrlbl'; endif; if cfact; label = ones(1,ni+1).*.label; endif; table = (blank~label)| (blank~b'); fmt = "*.*lf"~10~2; mask = zeros(1,cols(label)+1)| (0~ones(1,rows(b))); call printfm(table,mask,fmt); endif; preder = - c*b ; Erx = meanc(allr) ; if ergrf[1]; b2 = b; gosub grafit; /* graphs all. nonscaled only? */ endif; /* -------------------------------------------------------------------- */ /* 3. first stage standard errors, tests on b, 2nd-stage wghting matrix */ /* -------------------------------------------------------------------- */ j = 1; do while j <= rows(nwlags); s0 = finds0(f,nwlags[j]); W1 = inv(s0); varmat = 1/T*inv(C'W1*C); /* vcov of b */ se = diag(varmat)^(1/2); if prnt; format /RD 4,0; "Using" nwlags[j] "newey west lags "; table = ("s.e."~se')| ("t"~(b./se)'); mask = zeros(2,1)~ones(2,rows(se)); call printfm(table,mask,fmt); ""; /* chi 2 on coeffs */ gosub bchi2; /* "Predictions for each asset and std error in % "; let label = "er" "preder" "er-pred" "se"; $label'; 100*(erx~preder~(erx-preder)~diag(inv(W1)/T)^(1/2)); */ "rmse pricing error in percent, incl scaled "; ((meanc((preder-erx)^2))^(1/2))*100; "rmse pricing error, levels only "; ((meanc((preder[1:cols(rassetx)]- erx[1:cols(rassetx)])^2))^(1/2))*100; endif; j = j+1; endo; /* Hansen-Jaganathan Specification Error Test */ "";"Obsfact: Hansen-Jaganathan Spec Error Test; 1st stage"; format /RD 10,4; (lgt'inv((1/rows(allr))*allr'allr)*lgt)^.5; endif ; /* ------------------------- */ /* Part IV: Second stage GMM */ /* ------------------------- */ jtvec = zeros(rows(nwlags),1); dofvec = jtvec; if prnt; ""; "2. Second stage or iterated GMM "; endif; j = 1; do while j <= rows(nwlags); if prnt; " Using " nwlags[j] "newey west lags "; endif; format /RD 10,3; ""; /* provided w matrix section */ if (win /= 0) and (win /= 1); W1 = win; b1 = - inv(C'W1*C)*C'W1*Er; m1 = 1 + allfact*b1; f1 = (m1 .* allr); gT1 = meanc(f1); JT1 = T*gT1'W1*gT1; endif; /* iterated GMM */ if win == 1; lastw = eye(cols(f)); f1 = f; iter = 1; do while iter <= 100; /* -------------------------- */ /* 1. set up weighting matrix */ /* -------------------------- */ W1 = .5*lastw + .5*invpd(finds0(f1,nwlags)); lastw = W1; /* --------------------------------------------------- */ /* 2. with weighting matrix, set up JT & gT (from FOC) */ /* overid restrictions test */ /* --------------------------------------------------- */ b1 = - inv(C'W1*C)*C'W1*Er; m1 = 1 + allfact*b1; f1 = (m1 .* allr); gT1 = meanc(f1); JT1 = T*gT1'W1*gT1; if prnt; iter "th gmm iteration; JT this iter" Jt1; endif; /* convergence tests, break conditions */ /* test for .5% change in last 5 */ if iter == 1; oldchi = jt1; else; oldchi = jt1 | oldchi; endif; if iter > 5; test1 = abs(oldchi[1]-oldchi[5])/oldchi[1]; test2 = maxc(abs(oldchi[1:4]-oldchi[2:5])./oldchi[1]); if (test1 < 0.005) and (test2 < 0.005); if prnt; "ended at " iter " iterations" ; endif; break; endif; endif; if iter == 100; "obsfact: reached 100 gmm iterations without converging" ; "last % changes are: " test1*100 test2*100; endif; iter = iter + 1; endo; endif; /* two stage only */ /* -------------------------- */ /* 1. set up weighting matrix */ /* -------------------------- */ if win == 0; W1 = invpd(finds0(f,nwlags)); /* --------------------------------------------------- */ /* 2. with weighting matrix, set up JT & gT (from FOC) */ /* overid restrictions test */ /* --------------------------------------------------- */ b1 = - inv(C'W1*C)*C'W1*Er; m1 = 1 + allfact*b1; f1 = (m1 .* allr); gT1 = meanc(f1); JT1 = T*gT1'W1*gT1; endif; tneg_m = counts(m1,0)/rows(m1); jtvec[j] = jt1; dofvec[j] = cols(allr)-rows(b1); varmat = 1/T*inv(C'W1*C); se = diag(varmat)^(1/2); /* --------------------------------------------- */ /* 3. Print second stage estimates, if requested */ /* --------------------------------------------- */ if prnt; "Parameter estimates" ; /* note same varmat se as 1st */ let blank = " "; label = factlbl'; if xtrfact /= 0; label = label ~ xtrlbl'; endif; if cfact; label = ones(1,ni+1).*.label; endif; table = (blank ~ label)| /* stage GMM */ (" b:" ~ b1')| ("s.e.:" ~ se')| (" t:" ~ (b1./se)'); fmt = "*.*lf"~10~2; mask = zeros(1,cols(label)+1)| (zeros(3,1)~ones(3,rows(b1))); call printfm(table,mask,fmt); ""; b = b1; /* warning: b now useless. why is b,b1 distinct?*/ gosub bchi2; ""; /* ----------------------------------- */ /* 4. second stage overid restrictions */ /* ----------------------------------- */ "Overidentifying restrictions test" ; let label = "unrest."; chi2 = jt1; doff = cols(allr)-rows(b1); pval = 100*cdfchic(jt1,doff); if xtrfact /= 0; /* restricted system with this w matrix, no extra factor */ /* jtvec = obsfact(rassetx,factrs,xtrfact,inst, nwlags,ergrf,prnt,win,factlbl,xtrlbl,cfact,momcode,e, preder,erx,grphcode);*/ {jt,trsh,junk,pjunk,ejunk,gjunk} = obsfact(rassetx,factrs,0,inst, nwlags[j],0,0,w1,factlbl,xtrlbl,cfact,momcode,0,0,0,0, 2); /* note: it's pretty quick to do this by brute force instead */ label = label|"f only"; chi2 = chi2|jt; doff = doff|(cols(allr)-nf*(1+cfact*ni)); pval = pval|(100*(cdfchic(chi2[rows(chi2)],doff[rows(doff)]))); label = label|"f-unr."; if (jt-chi2[1])<0; chi2 = chi2|0; doff = doff|0; "Warning: bug, chi2 is going down when restricting"; else; chi2 = chi2|(jt-chi2[1]); doff = doff|(doff[rows(doff)]-doff[1]); endif; pval = pval|(100*(cdfchic(chi2[rows(chi2)],doff[rows(doff)]))); /* restricted system with this w matrix, only extra factors */ /* jtvec = obsfact(rassetx,factrs,xtrfact,inst, nwlags,ergrf,prnt,win,factlbl,xtrlbl,cfact,momcode, e,preder,erx,grphcode); */ {jt,trsh,junk,pjunk,ejunk,gjunk} = obsfact(rassetx,xtrfact,0,inst, nwlags[j],0,0,w1,xtrlbl,factlbl,cfact,momcode,0,0,0,0, 2); label = label|"xtronly"; chi2 = chi2|jt; doff = doff|(cols(allr)-nx*(cfact*ni+1)); pval = pval|(100*(cdfchic(chi2[rows(chi2)],doff[rows(doff)]))); label = label|"x-unr."; if (jt-chi2[1])<0; chi2 = chi2|0; doff = doff|0; else; chi2 = chi2|(jt-chi2[1]); doff = doff|(doff[rows(doff)]-doff[1]); endif; pval = pval|(100*(cdfchic(chi2[rows(chi2)],doff[rows(doff)]))); endif; format /RD 10,2; " " $label'; " JT" chi2'; format /RD 10,0; " DF" doff'; format /RD 10,3; " % p value" pval'; if ergrf[2]; preder = - c*b1 ; gosub grafit; endif; /* "Predictions for each asset and std error in % "; let label = "er" "preder" "er-pred" "se"; $label'; 100*(erx~preder~(erx-preder)~diag(inv(w)/T)^(1/2)); */ ""; endif; j = j+1; endo; /* Hansen-Jaganathan Specification Error Test */ if w_matrix == 0; gT1_=meanc((1+allfact*(-inv(C'inv((1/rows(allr))*allr'allr)*C)* C'inv((1/rows(allr))*allr'allr)*Er)).*allr); "";"Hansen-Jaganathan Spec Error Test -- 2nd stage"; " delta_W1 delta_W2"; ((gt1'inv((1/rows(allr))*allr'allr)*gt1)^.5)~((gt1_' inv((1/rows(allr))*allr'allr)*gt1_)^.5); elseif w_matrix == 1; /* w_matrix == 1 */ "";"Hansen-Jaganathan Spec Error Test -- 2nd stage"; " delta_W1"; ((gt1'inv((1/rows(allr))*allr'allr)*gt1)^.5); else; endif; if ergrf == 0; preder = 0; erx = 0; endif; retp(jtvec,dofvec,tneg_m,preder,erx,gt1); /* ---------------------------------------- */ /* Part V. Subroutines */ /* ---------------------------------------- */ /* ------------------------------------------------------------------ */ /* A. GRAFIT: sets up the publication graphics settings to construct */ /* E(r) vs. predicted E(r) graphs */ /* ------------------------------------------------------------------ */ grafit: graphset; _protate = 1; margin(0,0,0,3); _Plctrl = -1|0;_pltype = 6; ylabel("Mean Excess Return. % per qtr."); xlabel("Predicted Mean Excess Return, % per qtr."); /* xtics(0,3,.5,0); ytics(0,3,.5,0); */ if not ergrf[3]; _pscreen = 0; endif; if errbars; stderr = stdc(allr)/((rows(allr))^.5); _perrbar = 100*(preder*ones(1,3)~Erx~(Erx-stderr)~(Erx+stderr))~ ((6~2~1).*ones(rows(preder),3)); endif; _plegctl = 1; xtics(0,3.5,.5,1); ytics(0,3.5,.5,1); if grphcode == 3; preder2 = preder[1:cols(rassetx)]; preder = preder[cols(rassetx):rows(preder)]; erx2 = erx[1:cols(rassetx)]; erx = erx[cols(rassetx):rows(erx)]; endif; if (grphcode == 1) or (grphcode == 3) ; _Plctrl = -1|-1|0; _pstype = 8|3; _Psymsiz = 3|5; _plegstr = "Nonscaled Returns\000Scaled Returns\00045[o] Line"; xy(100*fill(preder2,fill(preder,(0|maxc(Erx)))), 100*fill(erx2,fill(erx,(0|maxc(erx))))); endif; if grphcode == 0; _pstype = 8; _plegstr = "Nonscaled Returns\00045[o] Line"; xy(100*fill(preder,(0|maxc(Erx))) , 100*fill(Erx,(0|maxc(Erx)))); endif; if grphcode == 3; preder = preder2|preder; erx = erx2|erx; endif; return; /* ------------------------------------------------------------------- */ /* B. BCHI2: Chi-Square tests of the parameter estimates. Testing */ /* whether (1) all factor estimates significant; (2) primary */ /* or additional factors are significant */ /* ------------------------------------------------------------------- */ bchi2: chi2 = b'inv(varmat)*b; doff = rows(b); pval = 100*cdfchic(chi2,doff); let label = "all b"; if nx > 0; /* compare model + extra */ select = seqb(1,nf,1); if cfact; k = 1; do while k <= cols(inst); select = select|seqb(k*(nf+nx)+1,k*(nf+nx)+nf,1); k = k+1; endo; endif; /* "model b select" select'; */ chi2 = chi2|b[select]'inv(varmat[select,select])*b[select]; doff = doff|rows(select); pval = pval | 100*cdfchic(chi2[rows(chi2)],doff[rows(doff)]); label = label|"model b"; select = seqb(nf+1,nf+nx,1); if cfact; k = 1; do while k <= cols(inst); select = select|seqb(k*(nf+nx)+nf+1,(k+1)*(nf+nx),1); k = k+1; endo; endif; /* "extra b select" select'; */ chi2 = chi2| b[select]'inv(varmat[select,select])*b[select]; doff = doff|rows(select); label = label|"extra b"; pval = pval | 100*cdfchic(chi2[rows(chi2)],doff[rows(doff)]); endif; if cfact and (nx == 0); /* compare model scaled + non */ select = seqb(1,nf,1); /* "nonscaled select" select'; */ chi2 = chi2| b[select]'inv(varmat[select,select])*b[select]; doff = doff|rows(select); label = label|"nonscaled"; pval = pval | 100*cdfchic(chi2[rows(chi2)],doff[rows(doff)]); select = seqb(nf+1,rows(b),1); /* "scaled select" select'; */ chi2 = chi2| b[select]'inv(varmat[select,select])*b[select]; doff = doff|rows(select); label = label|"scaled"; pval = pval | 100*cdfchic(chi2[rows(chi2)],doff[rows(doff)]); endif; "Chi-squared tests on b coefficients (%p is prob of higher chi-2)"; format /RD 10,4; label = blank | label; $label'; format /RD 10,2; " JT" chi2'; format /RD 10,0; " DF" doff'; format /RD 10,3; " % p value" pval'; return; endp;