Is effstim wrong when calculated in PHOTNU and PHOTLAM?

[pllim]

This does not make sense to me:

>>> import pysynphot as S
>>> S.__version__
'0.9.12'
>>> bp = S.ObsBandpass('acs,wfc1,f555w')
>>> sp = S.Vega
>>> obs = S.Observation(sp, bp)
>>> bp.pivot()
5360.938305731454
>>> obs.effstim()
1024.232011041325
>>> obs.effstim('photnu')
8.744735882086226e+28
>>> S.units.Photlam().ToPhotnu(bp.pivot(), obs.effstim())
9.818819146083892e-09
>>> obs.effstim('jy')
3630.073134871727
>>> S.units.Photlam().ToJy(bp.pivot(), obs.effstim())
3638.3434633950515

As shown above, effstim() calculated in Jy is not that far off compared to effstim() calculated in PHOTLAM and then converted to Jy around bandpass' pivot wavelength. However, that is not the case for PHOTNU. Why?

@deustua and @tddesjardins , since you have been looking into this from the INS side, what do you think? Is this a bug or am I not understanding something about effective stimulus in PHOTNU?

git blame shows that the Photnu.unitResponse() method used in the calculation in question was implemented by @vglaidler 11 years ago and has no changed since.