/* does means variances etc. on sp, crsp data */

library pgraph;
output file = prince.out reset;
graphs = 1;
convt = 1;

/* load data */

load x[] = \aaj\data\asia.dat;
x = reshape(x,rows(x)/35,35);
cal = x[.,1];
vwretd = x[.,2];
vwretx = x[.,4];

load x[] = \aaj\data\acti.dat;
x = reshape(x,rows(x)/21,21);
t90r = x[.,16];

/* 47 is obs 23 */
st = 23;
cal = cal[23:rows(cal)];
vwretd = vwretd[23:rows(vwretd)];
vwretx = vwretx[23:rows(vwretx)];
ustret = t90r[23:rows(t90r)];
T = rows(cal);
"Sample" cal[1] "to" cal[T];

/* 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;

"";
"one year returns on p/d, hh errors";
y = vwretd[2:T]-ustret[2:T];
x = ones(rows(y),1)~(1./dp[1:T-1]);
b = olsgmm(100*y,x,4);
"forecast" b'(1|(1./dp[T]));

"";
"two year returns on p/d, hh errors";
twoyrx =  100*((1+vwretd[2:T-1]).*(1+vwretd[3:T])-
             (1+ustret[2:T-1]).*(1+ustret[3:T]));
b2 = olsgmm(  twoyrx, ones(T-2,1)~1./dp[1:T-2],4);
"forecast" b2'(1|(1./dp[T]));

"";
"five year returns on p/d, hh errors";
fivxret = (1+vwretd[2:T-4]).*(1+vwretd[3:T-3]).*
             (1+vwretd[4:T-2]).*(1+vwretd[5:T-1]).*(1+vwretd[6:T])-
             (1+ustret[2:T-4]).*(1+ustret[3:T-3]).*
             (1+ustret[4:T-2]).*(1+ustret[5:T-1]).*(1+ustret[6:T]);
fivxret = 100*fivxret ;
b5 = olsgmm( fivxret, ones(T-5,1)~1./dp[1:T-5],10);
"forecast" b5'(1|(1./dp[T]));

graphjc;
title("Price/Dividend Ratio");
px = seqb(1947, 1998, 1);
py = 1./dp;
_pcolor = 15;
ytics(0,60,10,0);
xy(px,py);
graphps(px,py);
if convt;
 dos pqgrun graphic.tkf -c=1 -cf=princeton.eps;
endif;


/*
/* 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;