Overlaying/superimposing NBER recession dates on time series

I don't really know what you'd call it but I've always been curious how Econbrowser produces nice charts like this one on recession probabilities. I'm a SAS geek and this is how I coded it although I think that SAS should have an easier way of doing this without having to use an ANNOTATE data set which despite having used it quite a lot on project work is still a little bit of a mystery to me. Unfortunately, I don't know enough about HTML to render the code correctly so the input statement looks wrong.

data NBERdates;
   length PeakQ $ 20 TroughQ $ 20;
   informat Peak Trough ddmmyy10.;
   format Peak Trough date9.;
   input PeakQ 1-20 TroughQ 21-40 Peak Trough;
datalines;
February 1945(I)    October 1945 (IV)   1/2/1945        1/10/1945
November 1948(IV)   October 1949 (IV)   1/11/1948       1/10/1949
July 1953(II)       May 1954 (II)       1/7/1953        1/5/1954
August 1957(III)    April 1958 (II)     1/8/1957        1/4/1958
April 1960(II)      February 1961 (I)   1/4/1960        1/2/1961
December 1969(IV)   November 1970 (IV)  1/12/1969       1/11/1970
November 1973(IV)   March 1975 (I)      1/11/1973       1/3/1975
January 1980(I)     July 1980 (III)     1/1/1980        1/7/1980
July 1981(III)      November 1982 (IV)  1/7/1981        1/11/1982
July 1990(III)      March 1991(I)       1/7/1990        1/3/1991
March 2001(I)       November 2001 (IV)  1/3/2001        1/11/2001
;
run;

data NBERdates;
    set NBERdates;
    where Peak > '01Jan1947'd and Peak < '31Dec2000'd;
    PeakS = "'"put(Peak, date9.)"'d";
    TroughS = "'"put(Trough, date9.)"'d";
run;

filename gdp 'GDPC1.txt'; /* File downloaded from FRED St. Louis */
data GDP;
   infile gdp firstobs = 14;
   /* Read two variables: date and GDP */
   input date value;
   informat date yymmdd10.;
   format date yymmdd10.;
   lngdp = log(value);
   /* We will refer to Y instead of LNGDP for consistent notation */
   y=lngdp;
   /* Calculate the first difference of y */
   diffy = dif(y);
run;

%annomac;

data NBERDate1;
    set NBERdates;
    date = Peak; output;
    date = Trough; output;
run;

data anno;
    merge GDP(rename = (y=y1)) NBERdate1; by date;
    RETAIN YSYS XSYS '2';
    %bar(Peak,7.3,Trough,9.2,ligr,0,S);
run;

goptions ftext = swiss fontres = presentation htext = 0.9;
axis1 label = (angle=90);
axis2 order = ('1Jan1945'd to '30Jun2000'd by year5) minor = none;
symbol1 color = red interpol = j;
symbol2 color = black interpol = j;

proc gplot data = GDP;
plot lngdp * date /legend vaxis = axis1 haxis = axis2 annotate = anno chref=blue;
format date yyq6.;
run;
quit;