I'm trying to replicate a fan chart which shows a line along with +/- 1 and 2 standard errors. For context, this is an impulse response function. I've tried to replicate something similar (pictured below), but I'm still not getting a clean result. I'm open to any suggestions (I'm using Stata 15.1).

Here's my data.

Code:
* Example generated by -dataex-. To install: ssc install dataex
clear
input double hs float(period onese_plus onese_minus twose_plus twose_minus)
       0  1        0        0        0        0
 -.03425  2  -.02674  -.04176  -.01923  -.04927
-.044065  3 -.034615 -.053515 -.025165 -.062965
-.041269  4 -.030819 -.051719 -.020369 -.062169
  -.0446  5  -.03301  -.05619  -.02142  -.06778
-.042902  6 -.030302 -.055502 -.017702 -.068102
-.034948  7 -.021638 -.048258 -.008328 -.061568
-.029566  8 -.015646 -.043486 -.001726 -.057406
-.023311  9 -.009291 -.037331  .004729 -.051351
-.018024 10 -.003854 -.032194  .010316 -.046364
end
Here's what I've tried
Code:
twoway line hs period, color(red)   lwidth(thick) || ///
              rarea   onese_plus    onese_minus period,  fcolor(red%10)  lcolor(white)|| ///
              rarea   twose_plus     twose_minus period, fcolor(red%50)  lcolor(white) ///
                xtitle("Quarters") legend(region(lcolor(white)))
Here's a similar result I would like to achieve with the data I have.


Array