/*  part2.gpg
    does graph using standard errors from irstder.gpg

    index to graphs:
    fig1.tkf = VAR GNP, consumption            part2.gpg
    fig2.tkf = univariate GNP consumption
    fig3.tkf = spectral, GNP
    fig4.tkf = GNP ratios
    fig5.tkf = private GNP var                 part2.gpg
    fig6.tkf = private GNP univariate
    fig7.tkf = stock var                       part2.gpg
    fig8.tkf = stock univariate
    fig9.tkf = stock ratio
    fig10.tkf = bn trend
    fig11.tkf = bn big trend
    fig12.tkf = labor income var                part2.gpg
    fig13.tkf = labor income univariate
    fig14.tkf = labor income spectral

*/

library pgraph;
test = 0;         /* 1: look at all monte carlos to spot weird stuff */
mean = 0;         /* 1: center bars at m.c. mean; 0 center at xir */
                           /* try meadian, quantiles too */
errbars = 1;      /* put in error bars? */
dataset = 1;

if dataset == 1;
    load xirccv = xirccv1 ;
    load xircyv = xircyv1 ;
    load xirycv = xirycv1 ;
    load xiryyv = xiryyv1 ;
    load xir = xird1;
    horiz2 = 40;
endif;
if dataset == 2;
    load xirccv = xirccv2 ;
    load xircyv = xircyv2 ;
    load xirycv = xirycv2 ;
    load xiryyv = xiryyv2 ;
    load xir = xird2;
    horiz2 = 40;
endif;
if dataset == 3;
    load xirccv = xirccv3 ;
    load xircyv = xircyv3 ;
    load xirycv = xirycv3 ;
    load xiryyv = xiryyv3 ;
    load xir = xird3;
    horiz2 = 40;
endif;
if dataset == 4;
    load xirccv = xirccv4 ;
    load xircyv = xircyv4 ;
    load xirycv = xirycv4 ;
    load xiryyv = xiryyv4 ;
    load xir = xird4;
    horiz2 = 20;
endif;
if dataset == 5;
    load xirccv = xirccv5 ;
    load xircyv = xircyv5 ;
    load xirycv = xirycv5 ;
    load xiryyv = xiryyv5 ;
    load xir = xird5;
    horiz2 = 20;
endif;

select = 1|6|11|16|21|26|31|36;   /* keep conformable with irstder.gpg */

/* graph to see if ok */

if test;
    graphset;
    _Pltype = 6|(3*ones(rows(xirccv),1));
    _Pstype = 10;
    _Plctrl = 1|zeros(rows(xirccv),1);
    title("c,c shock");
    xy(select-1,xir[select,1]~xirccv');
    title("y,c shock");
    xy(select-1,xir[select,2]~xirycv');
    title("c,y shock");
    xy(select-1,xir[select,3]~xircyv');
    title("y,y shock");
    xy(select-1,xir[select,4]~xiryyv');
endif;

proc(1) = bigsort(x);             /* use to do means, quantiles */
    local i,xret;
    xret = sortc(x[.,1],1);
    if cols(x) > 1;
        i = 2;
        do while i <= cols(x);
            xret = xret~sortc(x[.,i],1);
            i = i+1;
        endo;
    endif;
    retp(xret);
endp;

proc(1) = bigindx(x,select);      /* matrix version of indexcat */
    local index,i;
    index = indexcat(x,select[1]);
    i = 2;
    do while i <= rows(select);
        index = index|indexcat(x,select[i]);
        i = i+1;
    endo;
    retp(index);
endp;


    graphset;
    _Pltype = 3|6|3|6;
    _Pstype = 3|3|2|2;
    _Plctrl = 5;
    fonts("simplex simgrma");
    ylabel("\201Response to 1 \202s\201 shock");
    xtics(0,horiz2,5,0);
    xlabel("Quarters after shock");
         xlbl = ""|"0"|
                ""|""|""|""|
                "10"|
                ""|""|""|""|
                "20"|
                ""|""|""|""|
                "30"|
                ""|""|""|""|
                "40";

    if dataset == 1;
        _Plegstr = "c, c shock\000y, c shock\000c, y shock\000y, y shock";
        _Plegctl=3|3|29|.9;
        _Ptek = "fig1.tkf";
         whichcc = 0|35;
         whichyc = 0|30;
         whichcy = 0|35;
         whichyy = 0|30;
         xoffset = 0|0;
         xtics(-2,horiz2,2,1);
         asclabel(xlbl,0);
         _Pmsgstr =
   "Fig. 1. Impulse-response function for consumption-GNP VAR (Table 1).  Lines"
$+
"\000give the response of consumption (c) and GNP (y) to one standard deviation"
$+
"\000shocks.  Bars show bootstrap one standard error bands.";
               /*  x y ht 0 0 10 */
    _Pmsgctl =  (1.7~6.7 ~.1~0~1~10)|
                (1.7~6.5 ~.1~0~1~10)|
                (1.7~6.3 ~.1~0~1~10);

    endif;

    if dataset == 2;
        @_Pline = 7~6~0~0~40~0~1~10; @       /* zero line */
        _Plegstr = "c, c shock\000e, c shock\000c, e shock\000e, e shock";
        _Plegctl=3|4|15|.8;
        ytics(-.25,1.25,.25,0);
        _ptek = "fig12.tkf";
         whichcc = 0|35;
         whichyc = 0|30;
         whichcy = 0|35;
         whichyy = 0|30;
         xoffset = -.75|0;
         xtics(-2,horiz2,2,1);
         asclabel(xlbl,0);
         _Pmsgstr =
"Fig. 12. Impulse-response function for consumption-labor income VAR (table 5)."
$+
"\000Lines give the response of consumption (c) and labor income (e) to one"
$+
" standard\000deviation "
$+
"shocks.  Bars show bootstrap one standard error bands.";

               /*  x y ht 0 0 10 */
    _Pmsgctl =  (1.7~6.7~.1~0~1~10)|
                (1.7~6.5 ~.1~0~1~10)|
                (1.7~6.3 ~.1~0~1~10);

    endif;

    if dataset == 3;
        _Plegstr = "c, c shock\000y, c shock\000c, y shock\000y, y shock";
        _Plegctl=3|4|25|.75;
        _Ptek = "fig5.tkf";
        ytics(-.5,1.5,.5,0);
         whichcc = 0|35;
         whichyc = 0|30;
         whichcy = 0|35;
         whichyy = 0|30;
         xoffset = 0|0;
         xtics(-2,horiz2,2,1);
         asclabel(xlbl,0);
         _Pmsgstr =
"Fig. 5. Impulse-response function for consumption-private GNP VAR (Table 3)."
$+
"\000Lines give the response of consumption (c) and private GNP (y) to one "
$+
"standard "
$+
"\000deviation shocks.  Bars show bootstrap one standard error bands."
;
               /*  x y ht 0 0 10 */
    _Pmsgctl =  (1.7~6.7~.1~0~1~10)|
                (1.7~6.5 ~.1~0~1~10)|
                (1.7~6.3 ~.1~0~1~10);
    endif;

    if (dataset == 4);
         xlabel("Years after shock");
        _Plegstr = "d, d shock\000p, d shock\000d, p shock\000p, p shock";
        _Ptek = "fig7.tkf";
        _Plegctl=3|3|14|4;
         ytics(-5,20,5,0);
         whichcc = 0|15;
         whichyc = 0|10;
         whichcy = 0|15;
         whichyy = 0|10;
         xoffset = .5|0;
         xtics(-1,horiz2,1,1);
         xlbl = ""|"0"|
                ""|""|""|""|
                "5"|
                ""|""|""|""|
                "10"|
                ""|""|""|""|
                "15"|
                ""|""|""|""|
                "20";
         asclabel(xlbl,0);
         _Pmsgstr =

"Fig. 7. Impulse-response function for dividend-price VAR (Table 4).  Lines "
$+
"\000give the response of dividends (d) and prices (p) to one standard"
$+
" deviation"
$+
"\000shocks.  Bars show bootstrap one standard error bands."
;
               /*  x y ht 0 0 10 */
    _Pmsgctl =  (1.7~6.7~.1~0~1~10)|
                (1.7~6.5 ~.1~0~1~10)|
                (1.7~6.3 ~.1~0~1~10);
    endif;

    if errbars;
        /* ERROR BARS */

        /* plot units: x left right y bot top linetype color.  */
        /* _Perrbar = (x ~ l ~ r  ~ y ~ b ~ t ~   lt   ~ 10);  */

        /* ----- */
        /*  c,c  */
        /* ----- */

        which = whichcc;
        index = bigindx(select,which+1);
        sdir = stdc(xirccv[.,index]);
        if mean;
            mir = meanc(xirccv[.,index]); /* to show biases in monte carlo */
        else;
            mir = xir[which+1,1];      /* to force through point estimate */
        endif;
       _Perrbar = (which+xoffset)~(which-.1+xoffset)~(which+.1+xoffset)~
                  mir~(mir-sdir)~(mir+sdir)~
                  (ones(rows(which),1)*(6~10));

        /* ----- */
        /*  y,c  */
        /* ----- */

        which = whichyc;         /* which irs to attach errors to, can be vec */
        index = bigindx(select,which+1);
        sdir = stdc(xirycv[.,index]);
        if mean;
            mir = meanc(xirycv[.,index]);       /*  to show biases */
        else;
            mir = xir[which+1,2];         /* to force through point estimate */
        endif;                                  /* try meadian, quantiles too */
       _Perrbar = _Perrbar|
                  (which-xoffset)~(which-.1-xoffset)~(which+.1-xoffset)~
                  mir~(mir-sdir)~(mir+sdir)~
                  (ones(rows(which),1)*(6~10));

        /* ----- */
        /*  c,y  */
        /* ----- */

        which = whichcy;         /* which irs to attach errors to, can be vec */
        index = bigindx(select,which+1);
        sdir = stdc(xircyv[.,index]);
        if mean;
            mir = meanc(xircyv[.,index]);
        else;
            mir = xir[which+1,3];
        endif;
       _Perrbar = _Perrbar|
                  which~(which-.1)~(which+.1)~
                  mir~(mir-sdir)~(mir+sdir)~
                  (ones(rows(which),1)*(6~10));

        /* ----- */
        /*  y,y  */
        /* ----- */

        which = whichyy;       /* which irs to attach errors to, can be vec */
        index = bigindx(select,which+1);
        sdir = stdc(xiryyv[.,index]);
        if mean;
            mir = meanc(xiryyv[.,index]);
        else;
            mir = xir[which+1,4];
        endif;
       _Perrbar = _Perrbar|
                  which~(which-.1)~(which+.1)~
                  mir~(mir-sdir)~(mir+sdir)~
                  (ones(rows(which),1)*(6~10));

    endif;

    xy(seqa(0,1,horiz2),xir[1:horiz2,.]);