/*              debtar1.prg
    takes maturity data from crsp.prg
    makes graphs that characterizes debt in "a cashless view"

    The back end
    estimates ar(1) surplus, debt process, for use in
   active debt policy exercise, but this was cut from the paper.
*/

/*
Do estimation of new bond sales. Note most sales are coupon bonds.

Play with model to reverse engineer desired properties

*/

library pgraph;
output file = debtar1.out reset;
graphs = 1;
rndseed 12345;
_olsres = 1;      /* tell ols to calculate residuals */
convt = 0;         /* create pic files */
gmtrc = 0;        /* fit geometric term structure or not */
bandw = 1;        /* 0: color 1: black and white */

/* --------- */
/* LOAD DATA */
/* --------- */

/* 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;
annual = 1;

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

if annual;
    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];
else;
    abase = abase[37:rows(abase)];    /* 1950:1 to 1009:12 */
endif;

/* load treasury data for comparison */

load x[] = fbfile.dat;
x = reshape(x,rows(x)/9,9);
x    = x[121:rows(x),.] ; /* keeps 68:1 on where all defined */
yrfb   = x[.,1];
mofb   = x[.,2];
FB   = x[.,3];  @ surplus           @
FBD  = x[.,4];  @ total debt        @
FBDP = x[.,5];  @ held by public    @
FBFP = x[.,6];  @ borrowing         @
FBO  = x[.,7];  @ outlays           @
FBOIN= x[.,8];  @ interest          @
FBR  = x[.,9];  @ receipts          @

sfb = chgfreq((fb+fboin)./punew[253:rows(punew)]*100,12,1,0); /* real surplus */
afb = chgfreq(fb./punew[253:rows(punew)]*100,12,1,0); /* conv. deficit*/
xfb = yrfb+mofb/12-1900;     /* NOTE 501231 is 51 here */
sfb = sfb[12:rows(sfb)];
xfb = xfb[12:rows(xfb)];
afb = afb[12:rows(afb)];

/* load Y,C data for macro graph */

load x[] = gdp.dat;
x = reshape(x,rows(x)/4,4);
yrq2 = x[.,1];
qq2 = x[.,2];
gd = x[.,3];
gdpq = x[.,4];

load x[] = cons.dat;
x = reshape(x,rows(x)/10,10);
yr3   = x[.,1];
qq3   = x[.,2];
GC   = x[.,3];
GCD  = x[.,4];
GCDQ = x[.,5];
GCN  = x[.,6];
GCNQ = x[.,7];
GCQ  = x[.,8];
GCS  = x[.,9];
GCSQ = x[.,10];

/* data starts in 59:1, but first two are zero so 59:3. Start it in 60 like
   the debt data */
xq = yrq2+qq2/4-1900;
cqa = chgfreq(gcnq+gcsq,4,1,0);  /* q obs of ann averages */
gdpqa = chgfreq(gdpq,4,1,0);
cqa = cqa[5:rows(cqa)];
gdpqa = gdpqa[5:rows(gdpqa)];    /* get rid of 59 */
xq = xq[5:rows(xq)];

cann = chgfreq(gcnq+gcsq,4,4,0)/4;
gdpann = chgfreq(gdpq,4,4,0)/4;
cann = cann[3:rows(cann)];
gdpann = gdpann[3:rows(gdpann)];

/* -- */
/* GO */
/* -- */

s = (-rev2-chabase)./apunew[2:rows(apunew)]*100;    /* real surplus */
                                /* NOTE price level is relative to 1, not 100.
                                   you MUST keep this convention below, use
                                   0.5 not 50 for price level */

 /* Chop all series up to 11 since few have private #*/

s = s[11:rows(s)];
matstr2 = matstr2[12:rows(matstr2),.];
abase = abase[12:rows(abase)];
matstr2[.,1] = matstr2[.,1]+abase;
mkt2 = mkt2[12:rows(mkt2)];

rev2 = rev2[11:rows(rev2)];
rev2c = rev2c[11:rows(rev2c)];
rev2d = rev2d[11:rows(rev2d)];
rev2m = rev2m[11:rows(rev2m)];
rev2n = rev2n[11:rows(rev2n)];
chabase = chabase[11:rows(chabase)];


/* eliminate maturities = columns with all zeros */
matstr2 = selif(matstr2',maxc(matstr2))';
ainfl = apunew[2:rows(apunew)]./apunew[1:rows(apunew)-1];
ainfl = ainfl[11:rows(ainfl)];
apunew = apunew[12:rows(apunew)];
T = rows(s);
dates = seqa(50+11,1,T);         /* NOTE 123150 is graphed as 50, not 51 */
                                 /* make sure higher frequency series are
                                    graphed conformably!!! */
/* load return to s calcluated s */

load sr = s;
sr = sr./apunew*100;

/* ------------------------ */
/* graph maturity structure */
/* ------------------------ */

if graphs ;


    /* s/c and infl */

    graphjc;
    if bandw;
       _pcolor = 15;
    endif;
    title("Surplus/Consumption and Inflation");
    px = dates;
    py = 100*(ainfl-1);
     py = (sr./cann/10)~ma(sr./cann/10,3)~py~ma(py,3);
                   /* ma makes centered 2 sided moving averages, see below */
    _pltype = 6|3|6|3;
    _plctrl = 1|0|0|0;

    _psymsiz = 2;
    _plegstr = "Surplus/consumption\000  moving average"$+
                "\000Inflation\000  moving average";
    _plegctl = 1|4|84|8;
    xtics(60,96,4,4);
    ytics(-6,15,3,3);
    ylabel("Percent");
    /*       1 typ xs ys xe ye 1=new co thick(0) */
    _pline = 1~6~61~0~97~0~1~15~0;
    if convt;
        _ptek = "sandpi.tkf";
    endif;
    xy(px,py);
    if convt;
        dos pqgrun sandpi.tkf -c=2 -cf=sandpi.pic;
    endif;

    proc(1) = ma(x,lag);
        local xbig;
        xbig = x*ones(1,2*lag+1);     /* copy x into many columns */
        xbig = xbig';                 /* data now in rows 1...T */
        xbig[1:lag,.] = shiftr(xbig[1:lag,.],rev(seqa(1,1,lag)),x[1]);
        xbig[lag+2:2*lag+1,.] =
                  shiftr(xbig[lag+2:2*lag+1,.],-seqa(1,1,lag),x[rows(x)]);
        xbig = meanc(xbig);
        retp(xbig);
   endp;



    /* s/c and y */

    graphjc;
    if bandw;
       _pcolor = 15;
    endif;
    title("Surplus and output");
    px = fill(dates,xq-1);       /* -1 so 96:4 lines up with 96 */
    py = gdpqa./cqa;
    py = 100*(py-meanc(py));
    py = fill(sr./cann/10,py);
    _pltype = 6|3;
    _plegstr = "Surplus/consumption\000Output/consumption - mean Y/C";
    _plegctl = 1|4|80|4;
    xtics(60,96,4,4);
    ytics(-8,8,4,4);
    ylabel("Percent");
    /*       1 typ xs ys xe ye 1=new co thick(0) */
    _pline = 1~6~61~0~97~0~1~15~0;
    if convt;
        _ptek = "sandcy.tkf";
    endif;
    xy(px,py);
    if convt;
        dos pqgrun sandcy.tkf -c=2 -cf=sandcy.pic;
    endif;



    /* nom, real debt growth and inflation */

    mkt = mkt2 + abase;
    nmkt = apunew.*mkt;
    dmkt = mkt[2:T]./mkt[1:T-1];
    dnmkt = nmkt[2:T]./nmkt[1:T-1];
    dsc = (sr[2:T]-sr[1:T-1])./cann[2:T]/1000;
    graphjc;
    fonts("simplex simgrma");
    begwind;
    if convt;
        _ptek = "debtgr.tkf";
    endif;
    _paxht = 0.25;
    _pnumht = 0.2;
    _ptitlht = 0.3;


    window(2,1,1);
    if bandw;
       _pcolor = 15;
    endif;
    px = dates[2:rows(dates)];
    py = 100*((dmkt~dnmkt)-1);
    xtics(60,96,4,4);
    ylabel("Percent");
    _Plegstr = "Real debt growth\000Nominal debt growth";
   _pltype = 6|3;
    _plegctl = 1|5|61|20;
    _pline = 1~6~61~0~97~0~1~15~0;
    ytics(-5,35,5,5);
    xy(px,py);

    nextwind;
    py = (100*(ainfl[2:rows(ainfl)]-1))~(100*sr[2:T]./cann[2:T]/1000);

   _pltype = 6|3;
    _psymsiz = 2;
    _plegstr = "Inflation\000Surplus/consumption";
    _plegctl = 1|5|61|8;
    /*       1 typ xs ys xe ye 1=new co thick(0) */
    _pline = 1~6~61~0~97~0~1~15~0;
    ytics(-10,15,5,5);
    xy(px,py);
    endwind;
    if convt;
        dos pqgrun debtgr.tkf -c=2 -cf=debtgr.pic;
    endif;

    /* print correlaiton statistics */

    "growth in debt value and d surplus / c ";
    "dates  real nominal inflation surplus";
    px~(100*((dmkt~dnmkt)-1))~(100*(ainfl[2:rows(ainfl)]-1))~(100*dsc);
    "";
    "correation of real debt growth, with change in s, s/c and s";
    sc = sr[2:T]./cann[2:T]/1000;

    corr(dmkt,dsc~sc~sr[2:rows(sr)]);

    /* graph history of surplus, conventional definition */

    graphjc;
    if bandw ;
        _pcolor = 15;
    endif;
    Title("Real Federal Surplus/Deficit");
    px = fill(dates,xfb-1);
    py = fill(sr/1000,afb);
    /*       1 typ xs ys xe ye 1=new co thick(0) */
    _pline = 1~6~61~0~97~0~1~15~0;
    ylabel("Billions 1987 $");
    xtics(60,97,5,5);
    ytics(-250,50,50,5);
    _plctrl = 1|0;
    @_plwidth = 3;@
    _pltype = 6|6;
    _pstype = 8|1;
    _psymsiz = 3;
    _plegstr =
              "Primary surplus, return calculation\000"
              $+ "Treasury surplus / deficit\000";
    _plegctl = 1|4|62|-240;
    if convt;
        _ptek = "surpls1.tkf";
        xy(px,py);
        dos pqgrun surpls1.tkf -cf=surpls1.pic -c=2;
    else;
        xy(px,py);
    endif;


    graphjc;
    if bandw ;
        _pcolor = 15;
    endif;
    Title("Three measures of primary surplus");
    px = fill((dates~dates),xfb-1);
                                   /* graph 123150 at 51 to be consistent */
                                   /* with monthly data */
    py = fill((sr~s)/1000,sfb);
    /*       1 typ xs ys xe ye 1=new co thick(0) */
    _pline = 1~6~61~0~97~0~1~15~0;
    ylabel("Billions 1987 $");
    xtics(60,97,5,5);
    ytics(-160,80,40,0);
    _plctrl = 1|1|0;
    @_plwidth = 3|1|1;@
    _pltype = 6|6|6;
    _pstype = 8|3|1;
    _psymsiz = 3;
    _plegstr =
              "Return calculation\000"
              $+ "Revenue calculation\000"
              $+"Net of interest surplus";
    _plegctl = 1|4|76|40;
    if convt;
        _ptek = "surplus.tkf";
        xy(px,py);
        dos pqgrun surplus.tkf -cf=surplus.pic -c=2;
    else;
        xy(px,py);
    endif;

    graphjc;
    if bandw ;
        _pcolor = 15;
    endif;
    title("Components of revenue from bond transactions");
    py = ((rev2n~rev2m~rev2c)./apunew*100)@~(-s)@~((rev2d~chabase)./apunew*100);
    py = py/1000;
    px = seqa(50+11,1,rows(py));
    _plegstr =
                   "Rev. = New bonds, bills"$+
               "\000     - Maturing bonds, bills"$+
               "\000     - Coupon payments"$+
               @"\000Total revenue = -s "$+ @
               "\000     + Ch. quantity outstanding"$+
               "\000     + Increase in base";
    _plegctl = 1|4|61|200;
    _psymsiz = 2;
    _plctrl = 1;
    ylabel("Billion 1987 $");
    xtics(60,96,4,4,);
    /*       1 typ xs ys xe ye 1=new co thick(0) */
    _pline = 1~6~61~0~97~0~1~15~0;
    if convt;
        _ptek = "compnts.tkf";
        xy(px,py);
        dos pqgrun compnts.tkf -cf=compnts.pic -c=2;
    else;
        xy(px,py);
    endif;


    graphjc;
    if bandw ;
        _pcolor = 15;
    endif;
    title("Maturity structure:  Zero-coupon face value ");
    px = dates;
    py = abase;
    py = py~(py[.,cols(py)]+matstr2[.,1]);
    py = py~(py[.,cols(py)]+sumc(matstr2[.,2:5]'));
    py = py~(py[.,cols(py)]+sumc(matstr2[.,6:10]'));
    py = py~(py[.,cols(py)]+sumc(matstr2[.,11:20]'));
    py = py~(py[.,cols(py)]+sumc(matstr2[.,21:cols(matstr2)]'));
    py = py./py[.,cols(py)]*100;
    _pmsgstr = "Monetary Base\0001 Year\0002-5 Year"$+
               "\0005-10 Years\00010-20 Years\00020+ years";
    filler = ones(6,1);
    _pmsgctl = (83|84|85|85|87|88)~
               (3 |25|55|75|85|95)~filler*(0.17~0~1~15~0);
@    _plegctl = 1|4|65|0.5;@
    _plctrl = 1;
    _psymsiz = 2;
    ylabel("Percent of total");
    xtics(60,96,4,4);
    if convt;
        _ptek = "matfrc.tkf";
        xy(px,py);
        dos pqgrun matfrc.tkf -cf=matfrc.pic -c=2;
    else;
        xy(px,py);
    endif;

    graphjc;
    if bandw ;
        _pcolor = 15;
    endif;
    title("Maturity structure, debt greater than 1 year maturity");
    px = seqa(0,1,cols(matstr2)-1);      /* maturities */
    py = ones(rows(px),1)*dates';
    pz = 100*(matstr2[.,2:cols(matstr2)]./sumc(matstr2'))';  /* don't graph 1 */
    viewxyz(80,-60,25);
    ztics(0,10,2,0);
    xlabel("Maturity");
    ylabel("Date"); zlabel("Percent");
    if convt;
        _ptek = "mat3d.tkf";
        xyz(px,py,pz);
        dos pqgrun mat3d.tkf -cf=mat3d.pic -c=2;
    else;
        xyz(px,py,pz);
    endif;

/*
do a 2 d version of this, like maturity structure

    graphjc;
    if bandw ;
        _pcolor = 15;
    endif;
    title("Sales/redemptions scaled by 0-1 year debt");
    px = seqa(1,1,cols(matstr2));
    py = ones(rows(px),1)*dates[2:rows(dates)]';
    pz = (matstr2[2:T,.]-(matstr2[1:T-1,2:cols(matstr2)]~zeros(T-1,1)))';
    pz = pz./matstr2[2:T,1]';
    viewxyz(40,-50,1);
    xlabel("maturity");
    ylabel("date");
    ztics(0,0.75,0.25,0);
    if convt;
        _ptek = "sale3d.tkf";
        xyz(px,py,pz);
        dos pqgrun sale3d.tkf -cf=sale3d.pic -c=2;
    else;
        xyz(px,py,pz);
    endif;
*/

endif;

/* characterize s for use in setting steady state */

"";
" s statistics"
"  mean" meanc(s);
"stddev" stdc(s);
"";

"";
"price statistics";
"average inflation (%)";
meandp = meanc(ln(apunew[2:T])-ln(apunew[1:T-1]));
100*meandp;

"";
"bond statistics";
"average growth in just maturing bonds (%)";
meandb1 = meanc(ln(matstr2[2:T,1]./matstr2[1:T-1,1]));
100*meandb1;
"";

/* --------------------------------------------------------  */
/* define steady state trends in bonds, price level, surplus */
/* --------------------------------------------------------  */

/* from here on, eliminate bonds more than 30 years */
matstr2 = matstr2[.,1:30];

rhv = ones(T,1)~seqa(1,1,T);
b = ln(matstr2[.,1])/rhv;
trdb0 = rhv*b;
Bbar = exp(trdb0[1]);
gb = exp(trdb0[2]-trdb0[1]);
B0bart = Bbar*gb^seqa(0,1,T);
"Bond steady state Bbar, gb " Bbar gb;

b = ln(apunew/100)/rhv;
trdp = rhv*b;
pbar = exp(trdp[1]);
gp = exp(trdp[2]-trdp[1]);
pbart = pbar*gp^seqa(0,1,T);
"Price steady state pbar, gp" pbar gp;

/* fit phi, steady state across bonds by eye; graph below */

relstr = (matstr2./matstr2[.,1]);
meanb = meanc(relstr);

if gmtrc;
    phi = 0.6;
else;
    phi = meanb;
endif;


r = 1.05;

/* infer s steady state from b,p */

gs = gb/gp;
if gmtrc;
    sbar = bbar/pbar*(r-gs)/(r-phi/gp);
else;
     sbar = bbar/pbar*(r-gs)/r*sumc(phi .* (1/r*gp)^seqa(0,1,rows(phi)));
endif;
sbart = sbar*gs^seqa(0,1,T);
"s steady state sbar gs " sbar gs ;

if graphs  ;

    graphjc;
    if bandw ;
        _pcolor = 15;
    endif;
    _plctrl = 1|0|0|0|0|0|0|0;
    _plegctl = 1;
    if gmtrc;
        _plegstr = "phi[j]\000grand mean\00060-75\00075-96";
        xy(seqa(0,1,cols(matstr2)),
              (phi^seqa(0,1,cols(matstr2)))~
               meanb~meanc(relstr[1:15,.])~meanc(relstr[16:36,.]));
    else;
        _plegstr = "phi]j[\000grand mean\00060-75\00075-96";
        xy(seqa(0,1,cols(matstr2)),
              phi~meanb~meanc(relstr[1:15,.])~meanc(relstr[16:36,.]));

        graphjc;
        if bandw ;
            _pcolor = 15;
        endif;
        _plctrl = 1;
        title("St.st. debt sales / Bbar gb^t; inferred from phi");
        xy(seqa(1,1,rows(phi)),phi*gb-(phi[2:rows(phi)]|0));

   endif;

    graphjc;
    title("surplus and steady state");
    _plctrl = 1|0;
    xy(dates,s~sbart);

    graphjc;
    title("price level and steady state");
    _plctrl = 1|0;
    xy(dates,(apunew/100)~pbart);

    graphjc;
    title("Just maturing bonds and steady state");
    _plctrl = 1|0;
    xy(dates,matstr2[.,1]~B0bart);

    graphjc;
    title("5,10,20 year bonds and steady state");
    _plctrl = 0;
    _pcolor = 1|2|3|1|2|3;
    if gmtrc;
        xy(dates,matstr2[.,6|11|21]~(B0bart*phi^(5~10~20)));
    else;
        xy(dates,matstr2[.,6|11|21]~(B0bart*phi[6|11|21]'));
    endif;



    graphjc;
    title("Bond sales and steady state");
    px = dates[2:T];
    py = b0bart[2:T].*(gb*phi'-(phi[2:rows(phi)]|0)'); /* st. st. sales */
    py2 = matstr2[2:T,.]-(matstr2[1:T-1,2:cols(matstr2)]~zeros(T-1,1));
    j = 1;
    do while j <= cols(py)/5;
        title(ftos((j-1)*5+1,"%*.*lf",3,0) $+ ftos(5*j,"%*.*lf",3,0));
        xy(px,py[.,(j-1)*5+1:j*5]~py2[.,(j-1)*5+1:j*5]);
        j = j+1;
    endo;




endif;


/* ------------------------------------------------- */
/* define deviations from steady state as in program */
/* ------------------------------------------------- */

ds = (s-sbart)./sbart;
dp = 100*pbart./apunew-1 ;            /* = 1/p - 1/pbar % 1/pbar  */
if gmtrc;
    db = matstr2./(b0bart*phi^seqa(0,1,cols(matstr2))')-1;
else;
    db = matstr2./(b0bart*phi')-1;
endif;


/* ---------------------------------------- */
/* run p ols on B,s to see if there is hope */
/* this produces a near perfect fit, as it should
   with N = T, but coefficients aren't unreasonable
   Compare with model predictions at some point
*/
/* ------------------------------------------------- */
/*
let __altnam = "const" "trend" "s"
       "B1" "B2" "B5"
       "B10" "B15" "B20"
       "-ln p";
__altnam = 0;
rhv = seqa(1,1,rows(apunew))~s~ln(matstr2[.,1:20]);
"";
"OLS regression of ln price level on bonds, s, to see if there is hope";
"";
{ vnam,m,b,stb,vc,stderr,sigma,cx,rsq,resid,dwstat } =
 ols("",-ln(apunew),rhv);

if graphs;
    graphjc;
    title("actual price level and fitted from bonds, surplus");
    xy(seqa(51+11,1,T),ln(apunew)~(-(ones(rows(rhv),1)~rhv)*b));
endif;
*/
/* ------------------------------------------------- */
/* Estimate joint surplus, debt rule                 */
/* since 35 years and 40 maturities, can't just run!
   instead, choose which coefficients carefully      */
/* ------------------------------------------------- */

/* pure variables */

cash = matstr2[.,1];
b1to5 = sumc(matstr2[.,2:6]');
blong= sumc(matstr2[.,7:cols(matstr2)]');
trends = ones(T-1,1)~seqa(1,1,T-1);
rhvcom = trends~s[1:T-1]~cash[1:T-1]~b1to5[1:T-1]~blong[1:T-1];
rhvcoms = trends~s[2:T]~s[1:T-1]~cash[1:T-1]~b1to5[1:T-1]~blong[1:T-1];
Nb = cols(matstr2);

/* deviations from "steady states" defined above */
dcash = db[.,1];
db1to5 = sumc(db[.,2:6]');
dblong= sumc(db[.,7:Nb]');
drhvcom = ds[1:T-1]~dcash[1:T-1]~db1to5[1:T-1]~dblong[1:T-1];
drhvcoms = ds[2:T]~ds[1:T-1]~dcash[1:T-1]~db1to5[1:T-1]~dblong[1:T-1];


/* --------------------------------------- */
/* use ols to look at elements of the VAR */
/* -------------------------------------- */

"s AR(1)";
let __altnam = "const" "s-1" "s";
call ols("",s[2:rows(s)],s[1:rows(s)-1]);
"";

"s AR(1) -- deviations";
let __altnam = "const" "s-1" "s";
call ols("",ds[2:T],ds[1:T-1]);
"";

"s AR(1) -- deviations, no constant";
__con = 0;
let __altnam =  "s-1" "s";
call ols("",ds[2:T],ds[1:T-1]);
"";
__con = 1;


"time trend";
let __altnam = "const" "t" "s";
{ vnam,m,b,stb,vc,stderr,sigma,cx,rsq,resid,dwstat } =
       ols("",s,seqa(1,1,rows(s)));

"";
"s  AR(1) with time trend";
let __altnam = "const" "t" "s-1" "s";
{ vnam,m,b,stb,vc,stderr,sigma,cx,rsq,resid,dwstat } =
       ols("",s[2:rows(s)],seqa(1,1,rows(s)-1)~s[1:rows(s)-1]);
"";
"s  AR(1) with time trend--deviations";
let __altnam = "const" "t" "s-1" "s";
{ vnam,m,b,stb,vc,stderr,sigma,cx,rsq,resid,dwstat } =
       ols("",ds[2:rows(s)],seqa(1,1,rows(s)-1)~ds[1:rows(s)-1]);

"";
" **************** S *************** ";
"";

let __altnam = "const" "trend" "s-1" "b0" "b1-5" "b6+" "s";
__output = 1;
"S on trend, lagged s, bonds ";
call ols("",s[2:T],rhvcom);

"";
let __altnam = "const" "s-1" "b0" "b1-5" "b6+" "s";
"Same, deviations S on trend, lagged s, bonds ";
call ols("",ds[2:T],drhvcom);

"";
"**************** B0 ******************* ";
"";

let __altnam = "const" "trend" "s" "s-1" "b0" "b1-5" "b6+" "b0";
__output = 1;
"";
"b0 on trend, s, lagged s, bonds ";
call ols("",cash[2:T],rhvcoms);

"";
let __altnam = "const" "s" "s-1" "b0" "b1-5" "b6+" "b0";
"Same deviaitons; b0 on trend, s, lagged s, bonds ";
call ols("",dcash[2:T],drhvcoms);


"";
"*************** Bshort ***************** ";
"";

let __altnam = "const" "trend" "s" "s-1" "b0" "b1-5" "b6+" "b1-5";
__output = 1;
"";
"b1-5 on trend, s, lagged s, bonds ";
call ols("",b1to5[2:T],rhvcoms);

"";
let __altnam = "const" "s" "s-1" "b0" "b1-5" "b6+" "b1-5";
"Same, deviations, b1-5 on trend, s, lagged s, bonds ";
call ols("",db1to5[2:T],drhvcoms);


"";
"****************** B long ********************* ";
"";
let __altnam = "const" "trend" "s" "s-1" "b0" "b1-5" "b6+" "blong";
"";
"blong on trend, s, lagged s, bonds ";
call ols("",blong[2:T],rhvcoms);
"";
let __altnam = "const" "s" "s-1" "b0" "b1-5" "b6+" "blong";
"Same, deviations: blong on trend, s, lagged s, bonds ";
call ols("",dblong[2:T],drhvcoms);


/* ------------ */
/* estimate VAR */
/* ------------ */


bigb = zeros(2+Nb,1+Nb);                /* place for coeffs in VAR format*/
bigb1 = zeros(cols(drhvcom)+3,1+Nb);    /* place for simple coeffs */

/* order: const trend s b------b */
/* indices of where to put things */
trendv = 1;
sv     = 2;
cashv  = 3;
b1to5v = 4|5|6|7|8;
nlong = cols(matstr2)-6;
blongv = seqa(9,1,nlong);

drhvcom = ones(rows(drhvcom),1)~drhvcom; /* add constant for / regressions */


bs = ds[2:T]/(drhvcom);
bigb1[1:cols(drhvcom),1] = bs;         /* no own, lagged coeff, obviously */
bigb[trendv,1] = bs[1];
bigb[sv,1] = bs[2];
bigb[cashv,1] = bs[3];
bigb[b1to5v,1] = bs[4]/5*ones(5,1);          /* spread sum coeffs evenly */
bigb[blongv,1] = bs[5]/nlong*ones(nlong,1);

let lbl = "    " "const" "s" "B0" "B1-5" "B5+";
$lbl';
" s " ;; bs';


lbl2 = lbl | "Bj" | "B1+1";
$lbl2';

j = 1;
do while j <= cols(matstr2)-1;
    if j > 1;
        bj = db[2:T,j]/(drhvcom~db[1:T-1,j:j+1]);
        bigb[2+j:2+j+1,j+1] = bj[6:7];          /* loading on own,next longer */
        bigb1[.,j+1] = bj;
    else;
        bj = db[2:T,j]/(drhvcom~db[1:T-1,j+1]);
        bigb[2+j,j+1] = bj[6];                    /* loading on own */
        bigb1[1:cols(drhvcom),j+1] = bj[1:rows(bj)-1];       /* no own,next */
        bigb1[cols(drhvcom)+2,j+1] = bj[rows(bj)];          /* j+1 goes here */
    endif;
   format /RD 6,2;
    "j ";; j ;; " bj ";; bj';

    bigb[trendv,j+1] = bigb[trendv,j+1]+bj[1];
    bigb[sv,j+1]     = bigb[sv,j+1]+bj[2];
    bigb[cashv,j+1]  = bigb[cashv,j+1]+bj[3];
    bigb[b1to5v,j+1] = bigb[b1to5v,j+1]+bj[4]/5*ones(5,1);
                             /* spread sum coeffs evenly */
    bigb[blongv,j+1] = bigb[blongv,j+1]+
                         bj[5]/cols(blong)*ones(cols(blong),1);
    j = j+1;
endo;

/* plot coefficients to develop smooth model of coeffs as fn of maturity */
/*

graphjc;
title("coefficients on s");
xlabel("maturity");
xy(seqa(1,1,cols(bigb1)-1),bigb1[2,2:cols(bigb1)]');
title("coefficients on B0");
xy(seqa(1,1,cols(bigb1)-1),bigb1[3,2:cols(bigb1)]');
title("coefficients on B1-5");
xy(seqa(1,1,cols(bigb1)-1),bigb1[4,2:cols(bigb1)]');
title("coefficients on B5+");
xy(seqa(1,1,cols(bigb1)-1),bigb1[5,2:cols(bigb1)]');

*/

output off;