/* does means variances etc. on sp, crsp data */ library pgraph; output file = meanvar.out reset; graphs = 1; convt = 0; domc = 0; /* se of mean monte carlo */ /* load data */ load x[] = asia.dat; x = reshape(x,rows(x)/35,35); cal = x[.,1]; vwretd = x[.,2]; vwretx = x[.,4]; ewretd = x[.,6]; ewretx = x[.,8]; decret = x[.,12|14|16|18|20|22|24|26|28|30]; load x[] = asbbi.dat; x = reshape(x,rows(x)/21,21); cstret = x[.,2]; gbtret = x[.,8]; cbtret = x[.,14]; ustret = x[.,16]; cpiret = x[.,18]; T = rows(x); load x[] = qsia.dat; /* 47:1 is obs 85 */ x = reshape(x,rows(x)/35,35); qcal = x[.,1]; qvwretd = x[.,2]; qvwretx = x[.,4]; qewretd = x[.,6]; qewretx = x[.,8]; qdecret = x[.,12|14|16|18|20|22|24|26|28|30]; load x[] = qsbbi.dat; x = reshape(x,rows(x)/21,21); qcstret = x[.,2]; qgbtret = x[.,8]; qcbtret = x[.,14]; qustret = x[.,16]; qcpiret = x[.,18]; qT = rows(x); load x[] = earndat.data; /* quarterly s&p data */ xe = reshape(x,rows(x)/5,5); ep = xe[.,2]./xe[.,5]; ep = ep[49:rows(ep)]; /* 49 = 4701 */ load x[] = inv1.dat; x = reshape(x,rows(x)/8,8); ifx = x[.,3]; /* total fixed */ in = x[.,5]; /* nonresidential */ ifx = ifx[3:rows(ifx)]; in = in[3:rows(in)]; /* 59:3 to 96:1 */ idat = x[.,1]+x[.,2]/4; idat = idat[3:rows(idat)]; load x[] = gdpq.dat; x = reshape(x,rows(x)/3,3); gdpq = x[3:rows(x),3]; /* clear x; */ /* means and variances */ basics = vwretd~cstret~ewretd~gbtret~cbtret~ustret; basics = (1+basics)./(1+cpiret)-1; let bsclbl = "VW" "SP" "EW" "GB" "CB" "US"; format /RD 7,2; "means, 1926-1996"; $bsclbl'; 100*meanc(basics)'; "means, 1947-1996"; $bsclbl'; 100*meanc(basics[22:T,.])'; "std dev" ; 100*stdc(basics[22:T,.])'; "std err sigma/root T"; 100*stdc(basics[22:T,.])'/(rows(basics[22:T,.])^0.5); "means, 1926-1996"; $bsclbl'; 100*meanc(basics)'; "std dev" ; 100*stdc(basics)'; "std err sigma/root T"; 100*stdc(basics)'/(rows(basics)^0.5); /* graph mean vs. sd. */ if graphs; graphjc; _pcolor = 15; _pcolor = 15; title("Mean vs. standard deviation of real returns, 1947-1996"); px1 = 100*stdc(basics[22:T,.]) ; px = fill(px1,100*stdc((1+decret[22:T,.])./(1+cpiret[22:T]))); py1 = 100*meanc(basics[22:T,.]); py = fill(py1,100*meanc((1+decret[22:T,.])./(1+cpiret[22:T])-1)); _pmsgstr = "Value weighted\000S&P 500\000Equally weighted" $+ "\000Govt bonds\000Corp bonds\000Treasury bill"; /* $+\000S1\000S2\000S3\000" $+"S4\000S5\000S6\000S7\000S8\000S9\000S10"; */ _pmsgctl = (px1[1:cols(basics),.]+(0.3))~ (py1[1:cols(basics),.]+(-.4|-.3|.2|-.2|0|0))~ (ones(cols(basics),1)*(0.12~0~1~15~0)); _plctrl = -1; _psymsiz = 3; _pstype = 10|2; xlabel("Standard deviation of return"); ylabel("Average return"); ""; "Mean vs. standard deviation of real returns"; format /RD 10,6; "x data"; px; "y data"; py; ""; xy(px,py); if convt; dos pqgrun graphic.tkf -c=2 -cf=msd.pic; endif; endif; /* table of std errors */ sig = 0.17; tim = 5|10|20|50; "table of standard errors as function of horizon"; tim~(100*(sig./tim^0.5))~(400*sig./tim^0.5); if domc; y = 100*(vwretd[22:T]-ustret[22:T]); yl = 100*(vwretd-ustret); i = 1; nt = 1000; do while i <= nt; ysim = y[trunc(rndu(rows(y),1)*rows(y))+1]; ysiml = yl[trunc(rndu(rows(yl),1)*rows(yl))+1]; if i == 1; means = meanc(ysim); meansl = meanc(ysiml); else; means = means|meanc(ysim); meansl = meansl|meanc(ysiml); endif; i = i+1; endo; "postwar" ; means = sortc(means,1); " mean y" meanc(y); "mean(mean(y))" meanc(means); "std (mean(y))" stdc(means); "quantiles 1 5 10 90 95 99" ; means[trunc((1|5|10|90|95|99)*nt/100+1)]'; "full sample"; meansl = sortc(meansl,1); " mean yl" meanc(yl); "mean(mean(yl))" meanc(meansl); "std (mean(yl))" stdc(meansl); "quantiles 1 5 10 90 95 99" ; meansl[trunc((1|5|10|90|95|99)*nt/100+1)]'; endif; /* returns on dp */ /* d/p = (R - Rx)/Rx = P(t)/P(t-1) + D(t)/P(t-1) - P(t)/P(t-1) = R/Rx - 1 */ /* d(t)/d(t-1) = d(t)/p(t) / d(t-1)/p(t-1) * p(t)/p(t-1) */ dp = (1+vwretd)./(1+vwretx)-1; dd = dp[2:T]./dp[1:T-1].*(1+vwretx[2:T])-1; dd = -999|dd; qdp = (1+qvwretd[1:qT-3]).*(1+qvwretd[2:qT-2]) .*(1+qvwretd[3:qT-1]).*(1+qvwretd[4:qT])./ ((1+qvwretx[1:qT-3]).*(1+qvwretx[2:qT-2]) .*(1+qvwretx[3:qT-1]).*(1+qvwretx[4:qT])) - 1; qdp = 0|0|0|qdp; /* pad out to same length */ ""; "one year returns on p/d, hh errors"; y = vwretd[23:T]-ustret[23:T]; x = ones(rows(y),1)~(1./dp[22:T-1]); b = olsgmm(100*y,x,4); ""; "two year returns on p/d, hh errors"; twoyrx = 100*((1+vwretd[23:T-1]).*(1+vwretd[24:T])- (1+ustret[23:T-1]).*(1+ustret[24:T])); b2 = olsgmm( twoyrx, ones(T-23,1)~1./dp[22:T-2],4); ""; "three year returns on p/d, hh errors"; threyrx = 100*((1+vwretd[23+0:T-2]).*(1+vwretd[23+1:T-1]).* (1+vwretd[23+2:T])- (1+ustret[23+0:T-2]).*(1+ustret[23+1:T-1]).* (1+ustret[23+2:T])); b3 = olsgmm( threyrx,ones(T-24,1)~1./dp[22:T-3],6); ""; "five year returns on p/d, hh errors"; fivxret = (1+vwretd[23+0:T-4]).*(1+vwretd[23+1:T-3]).* (1+vwretd[23+2:T-2]).*(1+vwretd[23+3:T-1]).*(1+vwretd[23+4:T])- (1+ustret[23+0:T-4]).*(1+ustret[23+1:T-3]).* (1+ustret[23+2:T-2]).*(1+ustret[23+3:T-1]).*(1+ustret[23+4:T]); fivxret = 100*fivxret ; b5 = olsgmm( fivxret, ones(T-26,1)~1./dp[22:T-5],10); "one year returns^2 on p/d, hh errors"; x = ones(rows(y),1)~(1./dp[22:T-1]); b1a = olsgmm((100*y)^2,x,4); ""; "two year returns^2 on p/d, hh errors"; b2a = olsgmm( twoyrx^2, ones(T-23,1)~1./dp[22:T-2],4); ""; "three year returns^2 on p/d, hh errors"; b3a = olsgmm( threyrx^2,ones(T-24,1)~1./dp[22:T-3],6); ""; "five year returns^2 on p/d, hh errors"; b5a = olsgmm( fivxret^2, ones(T-26,1)~1./dp[22:T-5],10); /* implied sigma. E(y^2) = ba*x; E(y) = b*x => sigma(y) = (ba*x-(b*x)^2)^0.5 */ if graphs; pdsim = seqb(10,30,1); one = ones(rows(pdsim),1); er = (one~pdsim)*(b~b2~b3~b5); sr = ((one~pdsim)*(b1a~b2a~b3a~b5a) - er^2); sr = sr./abs(sr) .* sqrt(abs(sr)); graphjc; title("Mean and std. dev of stock returns as a function of p/d"); _plegstr = "E1\000E2\000E3\000E5\000S1\000S2\000S3\000S5"; _plegctl = 1; /* xy(pdsim,er~sr); */ endif; /* Long horizon sharpe ratio statistics */ ""; "Mean and sd of long horizon excess simple returns"; "Horizon mean sd mean/h sd/h^0.5 sharpe sharpe/h^0.5 "; 1 100*meanc(y) 100*stdc(y) 100*meanc(y) 100*stdc(y) meanc(y)/stdc(y) meanc(y)/stdc(y); 2 meanc(twoyrx) stdc(twoyrx) meanc(twoyrx)/2 stdc(twoyrx)/(2^0.5) meanc(twoyrx)/stdc(twoyrx) meanc(twoyrx)/stdc(twoyrx)/(2^0.5); 3 meanc(threyrx) stdc(threyrx) meanc(threyrx)/3 stdc(threyrx)/(3^0.5) meanc(threyrx)/stdc(threyrx) meanc(threyrx)/stdc(threyrx)/(3^0.5); 5 meanc(fivxret) stdc(fivxret) meanc(fivxret)/5 stdc(fivxret)/(5^0.5) meanc(fivxret)/stdc(fivxret) meanc(fivxret)/stdc(fivxret)/(5^0.5); sevxret = (1+vwretd[23+0:T-6]).*(1+vwretd[23+1:T-5]).* (1+vwretd[23+2:T-4]).*(1+vwretd[23+3:T-3]).*(1+vwretd[23+4:T-2]) .*(1+vwretd[23+5:T-1]).*(1+vwretd[23+6:T])- (1+ustret[23+0:T-6]).*(1+ustret[23+1:T-5]).* (1+ustret[23+2:T-4]).*(1+ustret[23+3:T-3]).*(1+ustret[23+4:T-2]) .*(1+ustret[23+5:T-1]).*(1+ustret[23+6:T]); sevxret = 100*sevxret; tenxret = (1+vwretd[23+0:T-9]).*(1+vwretd[23+1:T-8]).* (1+vwretd[23+2:T-7]).*(1+vwretd[23+3:T-6]).*(1+vwretd[23+4:T-5]) .*(1+vwretd[23+5:T-4]).*(1+vwretd[23+6:T-3]).* (1+vwretd[23+7:T-2]).*(1+vwretd[23+8:T-1]).*(1+vwretd[23+9:T]) - (1+ustret[23+0:T-9]).*(1+ustret[23+1:T-8]).* (1+ustret[23+2:T-7]).*(1+ustret[23+3:T-6]).*(1+ustret[23+4:T-5]) .*(1+ustret[23+5:T-4]).*(1+ustret[23+6:T-3]).* (1+ustret[23+7:T-2]).*(1+ustret[23+8:T-1]).*(1+ustret[23+9:T]); tenxret = 100*tenxret; 7 meanc(sevxret) stdc(sevxret) meanc(sevxret)/7 stdc(sevxret)/(7^0.5) meanc(sevxret)/stdc(sevxret) meanc(sevxret)/stdc(sevxret)/(7^0.5); 10 meanc(tenxret) stdc(tenxret) meanc(tenxret)/10 stdc(tenxret)/(10^0.5) meanc(tenxret)/stdc(tenxret) meanc(tenxret)/stdc(tenxret)/(10^0.5); /* full sample */ "Full sample" ; yl = 100*(vwretd-ustret); twoyrxl = 100*((1+vwretd[1:T-1]).*(1+vwretd[2:T])- (1+ustret[1:T-1]).*(1+ustret[2:T])); threyrxl = 100*((1+vwretd[1:T-2]).*(1+vwretd[2:T-1]).* (1+vwretd[3:T])- (1+ustret[1:T-2]).*(1+ustret[2:T-1]).* (1+ustret[3:T])); fivxretl = 100*((1+vwretd[1:T-4]).*(1+vwretd[2:T-3]).* (1+vwretd[1+2:T-2]).*(1+vwretd[1+3:T-1]).*(1+vwretd[1+4:T])- (1+ustret[1+0:T-4]).*(1+ustret[1+1:T-3]).* (1+ustret[1+2:T-2]).*(1+ustret[1+3:T-1]).*(1+ustret[1+4:T])); 1 meanc(yl) stdc(yl) meanc(yl)/stdc(yl); 2 meanc(twoyrxl) stdc(twoyrxl) meanc(twoyrxl)/2 stdc(twoyrxl)/2^0.5 meanc(twoyrxl)/stdc(twoyrxl) meanc(twoyrxl)/stdc(twoyrxl)/2^0.5; 3 meanc(threyrxl) stdc(threyrxl) meanc(threyrxl)/3 stdc(threyrxl)/3^0.5 meanc(threyrxl)/stdc(threyrxl) meanc(threyrxl)/stdc(threyrxl)/3^0.5; 5 meanc(fivxretl) stdc(fivxretl) meanc(fivxretl)/5 stdc(fivxretl)/5^0.5 meanc(fivxretl)/stdc(fivxretl) meanc(fivxretl)/stdc(fivxretl)/5^0.5; /* logs */ "logs, postwar"; ""; yn = 100*(ln(1+vwretd[23:T])-ln(1+ustret[23:T])); twoyrxn = 100*(ln((1+vwretd[23:T-1]).*(1+vwretd[24:T]))- ln((1+ustret[23:T-1]).*(1+ustret[24:T]))); threyrxn = 100*(ln((1+vwretd[23+0:T-2]).*(1+vwretd[23+1:T-1]).* (1+vwretd[23+2:T]))- ln((1+ustret[23+0:T-2]).*(1+ustret[23+1:T-1]).* (1+ustret[23+2:T]))); fivxretn = 100*(ln((1+vwretd[23+0:T-4]).*(1+vwretd[23+1:T-3]).* (1+vwretd[23+2:T-2]).*(1+vwretd[23+3:T-1]).*(1+vwretd[23+4:T]))- ln((1+ustret[23+0:T-4]).*(1+ustret[23+1:T-3]).* (1+ustret[23+2:T-2]).*(1+ustret[23+3:T-1]).*(1+ustret[23+4:T]))); 1 meanc(yn) stdc(yn) meanc(yn)/stdc(yn); 2 meanc(twoyrxn) stdc(twoyrxn) meanc(twoyrxn)/2 stdc(twoyrxn)/2^0.5 meanc(twoyrxn)/stdc(twoyrxn) meanc(twoyrxn)/stdc(twoyrxn)/2^0.5 ; 3 meanc(threyrxn) stdc(threyrxn) meanc(threyrxn)/3 stdc(threyrxn)/3^0.5 meanc(threyrxn)/stdc(threyrxn) meanc(threyrxn)/stdc(threyrxn)/3^0.5; 5 meanc(fivxretn) stdc(fivxretn) meanc(fivxretn)/5 stdc(fivxretn)/5^0.5 meanc(fivxret)/stdc(fivxret) meanc(fivxret)/stdc(fivxret)/5^0.5; " full sample, not excess " ; "Full sample, not excess" ; yc = 100*(vwretd); twoyrxc = 100*((1+vwretd[1:T-1]).*(1+vwretd[2:T])-1); threyrxc = 100*((1+vwretd[1:T-2]).*(1+vwretd[2:T-1]).*(1+vwretd[3:T])-1); fivxretc = 100*((1+vwretd[1:T-4]).*(1+vwretd[2:T-3]).* (1+vwretd[1+2:T-2]).*(1+vwretd[1+3:T-1]).*(1+vwretd[1+4:T])-1); 1 meanc(yc) stdc(yc) meanc(yc)/stdc(yc); 2 meanc(twoyrxc) stdc(twoyrxc) meanc(twoyrxc)/2 stdc(twoyrxc)/2^0.5 meanc(twoyrxc)/stdc(twoyrxc) meanc(twoyrxc)/stdc(twoyrxc)/2^0.5; 3 meanc(threyrxc) stdc(threyrxc) meanc(threyrxc)/3 stdc(threyrxc)/3^0.5 meanc(threyrxc)/stdc(threyrxc) meanc(threyrxc)/stdc(threyrxc)/3^0.5; 5 meanc(fivxretc) stdc(fivxretc) meanc(fivxretc)/5 stdc(fivxretc)/5^0.5 meanc(fivxretc)/stdc(fivxretc) meanc(fivxretc)/stdc(fivxretc)/5^0.5; /* dividend growth regressions */ ""; "one year dividend growth on p/d, hh errors"; yd = dd[23:T]./(1+cpiret[23:T])-1; x = ones(rows(y),1)~(1./dp[22:T-1]); b = olsgmm(100*yd,x,4); ""; "two year dividend growth on p/d, hh errors"; call olsgmm( 100*( (1+dd[23:T-1]).*(1+dd[24:T])./ ( (1+cpiret[23:T-1]).*(1+cpiret[24:T]) )-1), ones(T-23,1)~1./dp[22:T-2],4); ""; "three year dividend growth on p/d, hh errors"; call olsgmm( 100*( (1+dd[23+0:T-2]).*(1+dd[23+1:T-1]).* (1+dd[23+2:T])./ ( (1+cpiret[23+0:T-2]).*(1+cpiret[23+1:T-1]).* (1+cpiret[23+2:T]) ) -1), ones(T-24,1)~1./dp[22:T-3],6); ""; "five year dividend growth on p/d, hh errors"; fivxd = 100* ( (1+dd[23+0:T-4]).*(1+dd[23+1:T-3]).* (1+dd[23+2:T-2]).*(1+dd[23+3:T-1]).*(1+dd[23+4:T])./( (1+cpiret[23+0:T-4]).*(1+cpiret[23+1:T-3]).* (1+cpiret[23+2:T-2]).*(1+cpiret[23+3:T-1]).*(1+cpiret[23+4:T]))-1); b5d = olsgmm( fivxd, ones(T-26,1)~1./dp[22:T-5],10); ""; if graphs; graphjc; _pcolor = 15; title("Price ratios"); px = seqb(47,96,1); py = 1./dp[22:T]; px = fill(px,seqb(47.25,96.75,0.25)); py = fill(py, 2*1./ep); _plegstr = "VW P/D\000 2 x S&P500 P/E"; _plegctl = 1|5|1955|2; _psymsiz = 2; _pstype = 3; _plctrl = 1|0; _plwidth = 5|1; xtics(1945,2000,5,0); ytics(0,50,10,0); px = 1900+px; ""; "Price ratios"; format /RD 10,6; " x data"; px; " y data" py; ""; xy(px,py); if convt; dos pqgrun graphic.tkf -c=2 -cf=pdpe.pic; endif; endif; /* prices and investment */ /* let k' = (1-delt)k + i; k'/i' = (1-delt)k/i i/i' + i/i' k'/i' = ((1-delt)k/i + 1) i/i' steady state (1-(1-d)i/i')k/i = i/i' k'/i' = i/i' / (1 - (1-d) i/i' ) */ delt = 0.05; /* per quarter */ ifxg = (ifx[2:rows(ifx)]./ifx[1:rows(ifx)-1]); ii1 = 1/meanc(ifxg); ki = ii1 / (1-(1-delt)*ii1); k = ki*(ifx[1]~in[1]); i = 1; do while i <= rows(ifx)-1; k = k|((1-delt)*k[rows(k),.]+(ifx[i]~in[i])); i = i+1; endo; "std dev of inv growth" 400*stdc(ln(ifx[2:rows(ifx)]./ifx[1:rows(ifx)-1])); if graphs; graphjc; title("p/d and investment"); px = idat; /* py = 600*(ifx)./k[.,1]; px = px~idat; */ py = (200*(ifx)./gdpq); /* px = fill(px,seqb(1947,1996,1)+1); /* plot stock returns at end of year */ py = fill(py,1./dp[22:T]); */ px = fill(px,seqb(1959.25,1997,0.25)); py = fill(py,1./qdp[133:rows(qdp)]); /* px = fill(px,seqb(1959.25,1997,0.25)); py2 = cumsumc(ln(qvwretx+1)); py = fill(py,py2[133:rows(py2)]); */ py = (py-meanc(py)')./stdc(py)'+(2~0); _pnum = 1|0; _plegstr = "I/Y\000VW P/D"; _plegctl = 1; _plctrl = 4|0; _pstype = 3; _psymsiz = 2; _pcolor = 15; xtics(1960,2000,5,5); ""; " P/D and investment "; format /RD 10,6; " x data"; px; " y data" py; ""; xy(px,py); if convt; dos pqgrun graphic.tkf -CF=pandi.pic -C=2; endif; endif; if graphs ; /* this version is not currently used in favor of version below that adds out of sample forecasts */ /* graphjc; _pcolor = 15; title("Actual and forecast five year returns"); px = seqb(47+5,96+5,1); py = (ones(T-21,1)~(1./dp[22:T]))*b5; /* forecast */ px = fill(px,seqb(47+5,91+5,1)); py = fill(py,fivxret); /* actual */ px = px + 1900; _plegstr = "Forecast\000Actual"; _plegctl = 1|5|1960|100; _plctrl = 0|1; _psymsiz = 2; _pltype = 6|4; _plwidth = 1|4; xlabel("Date"); ylabel("5 year excess return (%)"); xtics(1950,2005,5,0); ""; "Actual and forecast five year returns"; format /RD 10,6; " x data"; px; " y data" py; ""; xy(px,py); if convt; dos pqgrun graphic.tkf -c=2 -cf=5yrf&a.pic; endif; */ endif; /* look more at annual return regression */ /* forecast vs. actual plots etc. */ /* note since mean of pd[2:T] is higher than mean of pd[1:T-1] it is necessary to demean the variables before you start. Otherwise forecststs do not conver ge to mean pd and r */ y = vwretd[23:T]-ustret[23:T]; ydem = y - meanc(y); x = (1./dp[22:T-1]); xdem = x - meanc(x); b = ydem/xdem; " r on pd " b'; /* form long-term forecasts of annual returns from VAR */ y2 = (1./dp[23:T]); y2dem = y2 - meanc(y2); b2 = y2dem/xdem; "p/d ar(1)" b2'; /* r = -b1'- r-1 pd -b2'- pd-1 */ b2 = 0.98; "using pd ar(1)" b2; A = (0 ~b)| (0~b2); fcst = (vwretd[T]-ustret[T]-meanc(y))~(1./dp[T]-meanc(x)); j = 1; do while j <= 50; /* set forecast length here */ fcst = fcst|((A*(fcst[rows(fcst),.]'))'); j = j+1; endo; "forecast r and pd"; seqa(0,1,rows(fcst))~(100*(fcst[.,1]+meanc(y)))~(fcst[.,2]+meanc(x)); /* form long-term forecasts of annual returns from VAR */ if graphs ; graphjc; _pcolor = 15; title("Actual and forecast one year excess returns"); px = seqb(48,97,1); py = 100*(meanc(y)+((xdem*b)| ((1./dp[T]-meanc(x))*b)) ); /* in sample forecast */ /* last part adds 96 obs */ px = fill(px,seqb(48,96,1)); py = fill(py,100*y); /* actual returns */ px = fill(px,seqb(97,97+rows(fcst)-2,1)); /* out of sample forecast */ py = fill(py,100*(fcst[2:rows(fcst),1]+meanc(y))); /* px = fill(px,48|(97+rows(fcst)-1)); /* mean line */ py = fill(py,100*meanc(y)); */ _plegstr = "In-sample forecast\000Actual returns\000Out-of-sample forecast"; _plegctl = 1; px = px+1900; xtics(1945,2020,5,0); ytics(-35,45,5,0); xlabel("Date"); ylabel("Percent excess return"); _psymsiz = 2; _plctrl = 0|1|1; _pltype = 6|4|6 ; _plwidth = 1|4|1; ""; "Actual and forecast one year excess returns"; format /RD 10,6; " x data"; px; " y data" py; ""; xy(px,py); if convt; dos pqgrun graphic.tkf -c=2 -cf=lngfcst.pic; endif; endif; rhv = ones(rows(y),1)~x; b = y/rhv; fcst = rhv*b; if graphs ; graphjc; _pcolor = 15; title("Scatterplot of excess return (t+1) on p/d(t)"); _plctrl = -1|0|0; _psymsiz = 2; _pstype = 3; xlabel("p/d"); ylabel("Excess return"); _pmsgstr = "94\00095\00096\00097?"; _pmsgctl = ((1./dp[T-3]+0.3)~(100*y[rows(y)-2])~.15~0~1~15~0)| ((1./dp[T-2]+0.3)~(100*y[rows(y)-1])~.15~0~1~15~0)| ((1./dp[T-1]+0.3)~(100*y[rows(y)])~.15~0~1~15~0)| ((1./dp[T]+0.3)~(100*(1~(1./dp[T]))*b)~.15~0~1~15~0); xtics(10,45,5,0); py = 100*(y~(rhv*b)@~(0*x)@); py = py|(((100*(1~(1./dp[T]))*b)*(1~1))@~0@); px = x|(1./dp[T]); ""; "Scatterplot of excess return(t+1) on p/d(t)"; format /RD 10,6; " x data"; px; " y data" py; ""; xy(px,py); if convt; dos pqgrun graphic.tkf -c=2 -cf=sctplt.pic; endif; endif; /* same long term forecast graph with 5 year returns */ fivxret = (1+vwretd[23+0:T-4]).*(1+vwretd[23+1:T-3]).* (1+vwretd[23+2:T-2]).*(1+vwretd[23+3:T-1]).*(1+vwretd[23+4:T])- (1+ustret[23+0:T-4]).*(1+ustret[23+1:T-3]).* (1+ustret[23+2:T-2]).*(1+ustret[23+3:T-1]).*(1+ustret[23+4:T]); y5 = fivxret; y5dem = y5 - meanc(y5); x5 = (1./dp[22:T-5]); x5dem = x5 - meanc(x5); x5xdem = 1./dp[T-4:T]-meanc(x5); /* the extra right hand observations */ b5 = y5dem/x5dem; " r5 on pd " b5'; /* r = -b1'- r-1 pd -b2'- pd-1 */ "using pd ar(1)" b2; A5 = (0 ~b5)| (0~b2); fcst5 = (y5dem[rows(y5dem)])~(x5dem[rows(x5dem)]); /* last actual */ fcst5 = fcst5| ((x5xdem*b5)~x5xdem); /* to end of sample */ j = 1; do while j <= 50; /* set forecast length here */ fcst5 = fcst5|((A5*(fcst5[rows(fcst5),.]'))'); j = j+1; endo; fcst5 = fcst5[2:rows(fcst5),.]; /* elim. last actual */ "forecast r5 and pd"; seqa(0,1,rows(fcst5))~(100*(fcst5[.,1]+meanc(y5)))~(fcst5[.,2]+meanc(x5)); /* form long-term forecasts of annual returns from VAR */ if graphs ; graphjc; _pcolor = 15; title("Actual and forecast five year excess returns"); px = seqb(47+5,96,1); py = 100*(meanc(y5)+x5dem*b5); /* in sample forecast */ size(px) size(py); px = fill(px,seqb(47+5,96,1)); py = fill(py,100*y5); /* actual returns */ size(px) size(py); px = fill(px,seqb(96,96+rows(fcst5),1)); /* out of sample forecast */ py = fill(py,100*((meanc(y5)+x5dem[rows(x5dem)]*b5))| /* last in smpl */ (100*(fcst5[.,1]+meanc(y5)))); size(px) size(py); _plegstr = "In-sample forecast\000Actual returns\000Out-of-sample forecast"; _plegctl = 1|5|1960|120; px = px+1900; xtics(1950,2020,5,0); ytics(-60,180,30,0); xlabel("Date"); ylabel("Percent excess return"); _psymsiz = 2; _plctrl = 0|1|1; _pltype = 6|4|6 ; _plwidth = 1|4|1; ""; "Actual and forecast 5 year excess returns"; format /RD 10,6; " x data"; px; " y data" py; ""; xy(px,py); if convt; dos pqgrun graphic.tkf -c=2 -cf=lng5fcst.pic; endif; endif; proc(1) = olsgmm(y,x,lags); /* implements hh/newey west standard errors for ols regressions */ /* vcv = 1/T E(x(t)x(t)') ^-1 sum( w(j) E x(t)u(t)x(t-j)'u(t-j) ) E(x(t)x(t)') */ local T,b,e,xx1,s,j,vcv,stderr,s2,r2,olsse; T = rows(y); b = y/x; e = y-x*b; xx1 = inv(x'x/T); s = (e.*x)'(e.*x)/T; j = 1; do while j <= lags; s = s + (lags-j)/lags*( (e[1+j:T].*x[1+j:T,.])'(e[1:T-j].*x[1:T-j,.])/T + (e[1:T-j].*x[1:T-j,.])'(e[1+j:T].*x[1+j:T,.])/T ); j = j+1; endo; vcv = 1/T * xx1 * s * xx1; stderr = diag(vcv); stderr = stderr^0.5; s2 = e'e/T; olsse = s2*inv(x'x); olsse = diag(olsse); olsse = olsse^0.5; " coefficients" b'; " ols std errors" olsse'; "standard errors" stderr'; "t stat " b'./olsse'; " r2" 1-s2/stdc(y)^2; " resid correl" e[2:T]/e[1:T-1]; retp(b); endp;