Hi,
My colleague Alexey Tal (BSE expert) and I had a look at your input files. The parameters lead to GW and BSE results that are clearly underconverged with respect to the number of bands and k-points.
This is particularly the case for the BSE calculation since NBANDSV=1 and NBANDSO=1, while they should bet set to at least 12 in order to take into account all possible transitions up to 12 eV as in Fig. 3 of the paper (this can be deduced by looking at the GW quasiparticles energies). The number of k-points may also be too small.
Also, the GW band gap that we obtained is 7.1 eV, while 7.7. eV is reported in the paper. Thus, the number of k-points and NBANDS should be tested in the GW calculation.
Besides, it was also not possible to reproduce the spectrum that you generated. Did you really extract the imaginary part of the dielectric function printed in vasprun.xml or vaspout.h5? For future posts, please provide also the OUTCAR files of your calculations.