/* --------------------------------------------------------------------------
                           rettos.prg
   uses returns to find surplus. Idea is

   value(t)/value(t-1) = return + revenue/value(t-1)

it is easier and more accurate
to estimate the return based on the t-1 portfolio than to
properly figure out when bonds come due.

------------------------------------------------------------------------  */

library pgraph;
output file = rettos.out reset;
bandw = 1;          /* black and white?*/
convt = 0;         /* make pic files for import to scientific word?*/
testgrf = 0;      /* print out lots of diagnostic graphs */
                  /* 0 gives only graphs that are used in the paper */

/* ------------------------------------------------------------------- */
/* load series from bond analysis                                      */
/* ist obs of values is 501231, i.e. total for year 1950, date 1951    */
/* since rev is based on changes, its first year is 511231             */
/* ------------------------------------------------------------------- */

load rev1,rev2,face1,face2,mkt1,mkt2,matstr1,matstr2,facestr1,facestr2,
     rev2c,rev2d,rev2m,rev2n;

/* chop anything you're going to use to 601231 on (prices), 61 (returns)*/

T = rows(matstr2);
face2 = face2[11:T];
mkt2 = mkt2[11:T];
matstr2 = matstr2[11:T,.];
facestr2 = facestr2[11:T,.];

/* -------------------------------------------------------------- */
/* load fama french data, produce series of discount bond returns */
/* -------------------------------------------------------------- */

load x[] = famablis.dat;
x = reshape(x,rows(x)/6,6);
indx = seqa(7,12,rows(x)/12+1);
indx = selif(indx,indx .<= rows(x));   /* indices of end of year rows 1231xx */
x = x[indx,.];
fbdat = x[.,1];
p15 = x[.,2:6];                       /* bond prices 1-5 year maturity */
p15 = p15[9:rows(p15),.];             /* chop to 601231-now */


/* yields */
y15 = (100./p15)^(1./seqa(1,1,5)');

/* --------------------------------------------------------------------- */
/* impute p,returns for 6-30 year bonds by assuming same as 5 year yield */
/* --------------------------------------------------------------------- */

ybig = y15[.,5]*ones(1,cols(matstr2)-5);                  /* use 5 year yield */
                                                          /* for longer mat */
pbig = 100./(ybig^seqa(6,1,cols(matstr2)-5)');
pbig = p15~pbig;
ybig = y15~ybig;

/*  returns */
rbig = ((100*ones(rows(pbig)-1,1))~pbig[2:rows(pbig),1:cols(pbig)-1])
                  ./pbig[1:rows(pbig)-1,.];
                                         /* return = price now
                                          / price yesterday of one longer mat */



fbdat = seqa(60,1,rows(pbig));         /* dates for prices comparable with
                                          annual data, 123161 = 61 */

/* --------------------------------------- */
/* load 3 month t bill rate for short rate */
/* --------------------------------------- */

load x[] = fygn3.dat;
x = reshape(x,rows(x)/3,3);
fygn3 = 1+x[.,3]/100;
indx = seqa(12,12,rows(fygn3));
indx = selif(indx,indx.<=rows(fygn3));
y0 = fygn3[indx];                    /* annual , dec,dec, yield....*/
y0 = y0[14:rows(y0)];

indx = seqa(12,3,rows(fygn3));
indx = selif(indx,indx.<=rows(fygn3));
r0 = fygn3[indx];                    /* quarterly , dec,mar,jun,sep,dec, .*/
r0 = exp(chgfreq(ln(r0),4,4,0)/4);    /* annual returns */
r0 = r0[14:rows(r0)];

/* ------------------------------ */
/* load money base and price data */
/* ------------------------------ */

load xm[] = moneydat.dat;
xm  = reshape(xm,rows(xm)/7,7);
yr = xm[.,1];
Tm  = rows(yr);
mo = xm[.,2];
fmbase  = xm[.,5]*1000;                 /* crsp is millions, this is billions */
punew = xm[.,7];                        /* check that crsp is millions! */
clear xm;

abase = chgfreq(fmbase,1,12,0);    /* annual series, end of year */
apunew = chgfreq(punew,1,12,0);
abase = abase[4:rows(abase)];     /* 1950:12:1996:12 */
apunew = apunew[4:rows(apunew)];     /* 1950:12:1996:12 */
chabase = abase[2:rows(abase)]-abase[1:rows(abase)-1];

abase = abase[11:rows(abase)] ;
apunew = apunew[11:rows(apunew)];
chabase = chabase[10:rows(chabase)];    /* all now 611231-now */

infl = apunew[2:rows(apunew)]./apunew[1:rows(apunew)-1]; /* 62 on */
mkt = mkt2+abase;                  /* total market value privately held inc M*/


/* ------------------------------------ */
/* graphs to check that it's all sensible */
/* ------------------------------------ */

if testgrf;
    graphjc;
    title("year t-bill returns and FB returns");
    _plegstr = "0\0001\0002\0003\0004\0005\00010";
    _plegctl = 1;
    _plctrl = 1;
    xy(fbdat[2:rows(fbdat)],r0~rbig[.,1|2|3|4|5|10]);

    title("december t bill rate and FB yields");
    _plegstr = "0\0001\0002\0003\0004\0005";
    _plegctl = 1;
    _plctrl = 1;
    xy(fbdat,y0~ybig[.,1|2|3|4|5])


    title("FB prices");
    _plegstr = "1\0002\0003\0004\0005";
    _plegctl = 1;
    _plctrl = 1;
    xy(fbdat,pbig[.,1|2|3|4|5]);
endif;


/* ------------------------------------------ */
/* compate value at FB prices to exact value  */
/* ------------------------------------------ */

tvfb1 = sumc((matstr2.*pbig/100)'); /* total value at FB prices */
                                    /* 0-1 debt is 1 year, 1-2 is 2 year etc*/
fbprice = ones(rows(matstr2),1) ~ (pbig[.,1:cols(pbig)-1]/100);
tvfb2 = sumc((matstr2.*fbprice)');
                                    /* total value at FB prices */
                                    /* 0-1 debt is 0 year, 1-2 is 1 year etc*/
                                    /* this convention is closest to actual */

tvfb1 = tvfb1 + abase;              /* add base to total value */
tvfb2 = tvfb2 + abase;

if testgrf;
    graphjc;
    title("Comparing value at FB prices to actual; verifies FB is ok approx"$+
           "\LTotal nominal market value of debt");
    _plegstr = "FB, 0-1 = 1\000FB, 0-1 = 0\000actual";
    _plegctl = 1;
    py = tvfb1~tvfb2~mkt;
    xy(fbdat,py);
    title("Comparing value at FB prices to actual; verifies FB is ok approx"$+
           "\LGrowth in total nominal market value of debt");
    xy(fbdat[2:rows(fbdat)],py[2:rows(py),.]./py[1:rows(py)-1,.]);
endif;

/* -------------------------------------------------- */
/* compute rate of return on government bond portfolio */
/* return for year t (december, t) is
   sum return_i (t-1 ->t) * bonds(t-1,i) * price(t-1,i)/value(t-1)
   use FB value for consistency
   use T bill return for 0-1 years.                  */
/* -------------------------------------------------- */

T = rows(matstr2);
wts = (abase[1:T-1,.]~matstr2[1:T-1,.]).*             /* base ~ maturities */
       (ones(T-1,1)~fbprice[1:T-1,.])                    /* prices            */
             ./tvfb2[1:T-1];
rgov = (ones(rows(r0),1)~r0~rbig[.,1:cols(rbig)-1]).* wts;
rgov = sumc(rgov');


if testgrf;
    graphjc;
    title("govt bond return and actual, FB value growth");
 _plegstr = "govt bond return\000actual value growth\000FB prices value growth";
    _plegctl = 1;
    xy(fbdat[2:T],rgov~(tvfb2[2:T]./tvfb2[1:T-1])~
         (mkt[2:T]./mkt[1:T-1]));
endif;



graphjc;
if bandw;
    _pcolor = 15;
endif;
title("Nominal government bond returns");
_plegstr = "Gov't average\0003 Month\0005 Year";
_psymsiz = 3|3|1;
_plctrl = 1|0|0|0|0|0|0;
_plegctl = 1|4|80|-4;
_pltype = 6|6|3;
ytics(-5,20,5,0);
ylabel("Percent return");
xtics(60,96,4,4);
if convt;
    _ptek = "rongov.tkf";
endif;
xy(fbdat[2:T],100*((rgov~r0~rbig[.,5])-1));
if convt;
    dos pqgrun rongov.tkf -c=2 -cf=rongov.pic;
endif;



/* -------------------------------------------------- */
/* finally, compute surplus as residual,
   V(t)/V(t-1) = r(t-1 -> t) - s(t)/V(t-1)
   v(t-1)*r(t-1 -> t) - V(t) = s(t)                  */
/* -------------------------------------------------- */

s = (mkt[1:T-1]).*rgov - (mkt[2:T]);


graphjc;
title("Real Surplus");
_plegstr = "Inferred from value growth, return\000direct, estimated revenue";
_plegctl = 1|4|62|50;
_pltype = 6|1;
xy(fbdat[2:T],(s~(-rev2[11:rows(rev2)]-chabase[2:rows(chabase)]))./
               apunew[2:rows(apunew)]*100/1000);

graphjc;
title("Real surplus");
ylabel("$ Billion");
xy(fbdat[2:T],s./apunew[2:rows(apunew)]*100/1000);


graphjc;
if bandw;
    _pcolor = 15;
endif;
title("Real government bond return, surplus and debt growth");
_plegstr = "bond return\000debt growth\000surplus/debt";
_plegctl = 1;

/* definition implies
s./mkt[1:T-1]./infl = rgov./infl - (mkt[2:T]./mkt[1:T-1])./infl;
*/

py = (rgov./infl)~((mkt[2:T]./mkt[1:T-1])./infl)~(s./mkt[1:T-1]./infl);
py = ((py[.,1:2]-1)~py[.,3])*100;
"Correlations of return, debt growth, surplus/debt";
corr(py,py);
_psymsiz = 2;
_plctrl = 0|0|1;
_pltype = 6|3;
ylabel("percent");
    /*       1 typ xs ys xe ye 1=new co thick(0) */
    _pline = 1~6~62~0~96~0~1~15~0;
xtics(60,96,4,4);
if convt;
    _ptek = "randebt.tkf";
endif;
xy(fbdat[2:T],py);
if convt;
    dos pqgrun randebt.tkf -c=2 -cf=randebt.pic;
endif;


save s,rgov;

output off;

output file = aggreg2.txt reset;
"Debt aggregates";
"Base = monetary base (citibase)";
"Market value = privately held debt + monetary base";
"Rate of return is weighted average of fama-bliss zero coupon returns "
" and three month rate, weights = portfolio shares at beginning of year,";
" Units = percent";
"Surplus[t] = value[t-1]*rate of return[t] - value[t]";
"CPI = Dec CPI (citibase)";
"";
"      date        base        value          return     surplus         CPI";
outwidth 250;
format /RD 12,4;
seqa(60,1,rows(mkt))~abase~mkt~(-99|(100*(rgov-1)))~(-99|s)~apunew;


output off;