/* part3.gpg
    univariate impulse response with standard errors

    index to graphs:

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

*/

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 = 4;

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;

if dataset == 1;
    load xir2v = xir2v1 ;
    load xir1v = xir1v1 ;
    load xir1 = xir1d1;
    load xir2 = xir2d1;
    horiz2 = 40;
endif;
if dataset == 2;
    load xir2v = xir2v2 ;
    load xir1v = xir1v2 ;
    load xir1 = xir1d2;
    load xir2 = xir2d2;
    horiz2 = 40;
endif;
if dataset == 3;
    load xir2v = xir2v3 ;
    load xir1v = xir1v3 ;
    load xir1 = xir1d3;
    load xir2 = xir2d3;
    horiz2 = 40;
endif;
if dataset == 4;
    load xir2v = xir2v4 ;
    load xir1v = xir1v4 ;
    load xir1 = xir1d4;
    load xir2 = xir2d4;
    horiz2 = 20;
endif;
if dataset == 5;
    load xir2v = xir2v5 ;
    load xir1v = xir1v5 ;
    load xir1 = xir1d5;
    load xir2 = xir2d5;
    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(xir1v),1));
    _Pstype = 10;
    _Plctrl = 1|zeros(rows(xir1v),1);
    title("univariate ar");
    xy(select-1,xir1[select]~xir1v');
    title("var");
    xy(select-1,xir2[select]~xir2v');
endif;

    graphset;
    _Pltype = 3|6;
    ylabel("Response to a unit shock");
    px = seqa(0,1,horiz2);
    py = xir1[1:horiz2]~xir2[1:horiz2];

    if (dataset == 1);
        ytics(0,1.8,.2,0);
        _Ptek = "fig2.tkf";
        xlabel("Quarters after a unit shock");
        xtics(0,horiz2,5,0);
        _Pmsgstr =
    "Fig. 2. Univariate GNP impulse-response functions.  The 'univariate' "
$+
"\000response is esimated from a regression of GNP growth on past GNP growth."
$+
"\000The 'VAR' response is the univariate GNP response implied by the"
$+
"\000consumption-GNP VAR.  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)|
                (1.7~6.1 ~.1~0~1~10);


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

        which = 20|30;    /* which irs to attach errors to, can be vec */
        index = bigindx(select,which+1);
        sdir = stdc(xir1v[.,index]);
        /* mir = meanc(xir1v[.,index]); */ /*to display biases in monte carlo */
        mir = xir1[which+1,1];             /* to force through point estimate */
                                           /* try meadian, quantiles too */
       _Perrbar = which~(which-.1)~(which+.1)~
                  mir~(mir-sdir)~(mir+sdir)~
                  (ones(rows(which),1)*(6~10));

        which = 25|35;    /* which irs to attach errors to, can be vec */
        index = bigindx(select,which+1);
        sdir = stdc(xir2v[.,index]);
        /* mir = meanc(xir2v[.,index]); */ /*to display biases in monte carlo */
        mir = xir2[which+1,1];             /* to force through point estimate */
                                           /* try meadian, quantiles too */
       _Perrbar = _Perrbar |
                  which~(which-.1)~(which+.1)~
                  mir~(mir-sdir)~(mir+sdir)~
                  (ones(rows(which),1)*(6~10));


    endif;



    if dataset == 2;
        xlabel("Quarters after a unit shock");
        ytics(0,1.4,.2,0);
        _Ptek = "fig13.tkf";
        which = 20|30;    /* which irs to attach errors to, can be vec */
        index = bigindx(select,which+1);
        sdir = stdc(xir1v[.,index]);
        /* mir = meanc(xir1v[.,index]); */ /*to display biases in monte carlo */
        mir = xir1[which+1,1];             /* to force through point estimate */
                                           /* try meadian, quantiles too */
       _Perrbar = which~(which-.1)~(which+.1)~
                  mir~(mir-sdir)~(mir+sdir)~
                  (ones(rows(which),1)*(6~10));

        which = 25|35;    /* which irs to attach errors to, can be vec */
        index = bigindx(select,which+1);
        sdir = stdc(xir2v[.,index]);
        /* mir = meanc(xir2v[.,index]); */ /*to display biases in monte carlo */
        mir = xir2[which+1,1];             /* to force through point estimate */
                                           /* try meadian, quantiles too */
       _Perrbar = _Perrbar |
                  which~(which-.1)~(which+.1)~
                  mir~(mir-sdir)~(mir+sdir)~
                  (ones(rows(which),1)*(6~10));

_Pmsgstr =
"Fig. 13. Univariate labor income impulse-response functions.  The 'univariate'"
$+
"\000response is esimated from a regression of income growth on past income"
$+
" growth."
$+
"\000The 'VAR' response is the univariate response implied by the consumption -"
$+
"\000labor income VAR.  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)|
                (1.7~6.1 ~.1~0~1~10);


    endif;


    if dataset == 3;
        xlabel("Quarters after a unit shock");
        ytics(0,1.8,.2,0);
        _Ptek = "fig6.tkf";
        xtics(0,horiz2,5,0);
        _Pmsgstr =
 "Fig. 6. Univariate private GNP impulse-response functions.  The 'univariate'"
$+
"\000response is esimated from a regression of GNP growth on past GNP growth."
$+
"\000The 'VAR' response is the univariate GNP response implied by the"
$+
 " consumption - "
$+
"\000private GNP VAR.  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)|
                (1.7~6.1 ~.1~0~1~10);

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

        which = 20|30;    /* which irs to attach errors to, can be vec */
        index = bigindx(select,which+1);
        sdir = stdc(xir1v[.,index]);
        /* mir = meanc(xir1v[.,index]); */ /*to display biases in monte carlo */
        mir = xir1[which+1,1];             /* to force through point estimate */
                                           /* try meadian, quantiles too */
       _Perrbar = which~(which-.1)~(which+.1)~
                  mir~(mir-sdir)~(mir+sdir)~
                  (ones(rows(which),1)*(6~10));

        which = 25|35;    /* which irs to attach errors to, can be vec */
        index = bigindx(select,which+1);
        sdir = stdc(xir2v[.,index]);
        /* mir = meanc(xir2v[.,index]); */ /*to display biases in monte carlo */
        mir = xir2[which+1,1];             /* to force through point estimate */
                                           /* try meadian, quantiles too */
       _Perrbar = _Perrbar |
                  which~(which-.1)~(which+.1)~
                  mir~(mir-sdir)~(mir+sdir)~
                  (ones(rows(which),1)*(6~10));

    endif;
    if dataset == 4;
         xlabel("Years after a unit shock");
         ytics(0,1.25,.25,0);
        _Ptek = "fig8.tkf";
        _Pmsgstr =
 "Fig. 8. Univariate stock price impulse-response functions.  The 'univariate'"
$+
"\000response is estimated from a regression of returns on past returns."
$+
"\000The 'VAR' response is the univariate price response implied by the"
$+
" dividend - "
$+
"\000price VAR.  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)|
                (1.7~6.1 ~.1~0~1~10);

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

        which = 10;    /* which irs to attach errors to, can be vec */
        index = bigindx(select,which+1);
        sdir = stdc(xir1v[.,index]);
        /* mir = meanc(xir1v[.,index]); */ /*to display biases in monte carlo */
        mir = xir1[which+1,1];             /* to force through point estimate */
                                           /* try meadian, quantiles too */
       _Perrbar = which~(which-.1)~(which+.1)~
                  mir~(mir-sdir)~(mir+sdir)~
                  (ones(rows(which),1)*(6~10));

        which = 15;    /* which irs to attach errors to, can be vec */
        index = bigindx(select,which+1);
        sdir = stdc(xir2v[.,index]);
        /* mir = meanc(xir2v[.,index]); */ /*to display biases in monte carlo */
        mir = xir2[which+1,1];             /* to force through point estimate */
                                           /* try meadian, quantiles too */
       _Perrbar = _Perrbar |
                  which~(which-.1)~(which+.1)~
                  mir~(mir-sdir)~(mir+sdir)~
                  (ones(rows(which),1)*(6~10));

    endif;

    if dataset == 5;
        xlabel("Years after a unit shock");
        ytics(0,1.0,.25,0);
        _Ptek = "stk472.tkf";
    endif;

    _Plegctl = 1;
    _Plegstr = "Univariate\000VAR" ;
    @title("y impulse responses");@
    xy(px,py);