This has to do with the way the IRFs are computed for 2nd order. I turn to it below, but let me first confirm that the simulation is correct. With option periods=2100, you create one Monte Carlos simulation for 2100 periods.
I added variable EPV=PV(+1) to store the expected value of next period PV.
For the first three periods, we have indeed correct results (the two first periods are steady state values required by the initial lags on two periods). If you type
[PV(1:5) DIV(1:5) beta*EPV(1:5) DIV(1:5)+beta*EPV(1:5)]
after running the model, you get indeed
2.15132113253737 0.17210569060299 1.97921544193438 2.15132113253737
2.15132113253737 0.17210569060299 1.97921544193438 2.15132113253737
1.95654202503276 0.12947751161443 1.82706451341762 1.95654202503206
2.11397231142171 0.16609479510913 1.94787751631264 2.11397231142178
2.06428867370417 0.15028749532325 1.91400117838059 2.06428867370384
The problem with the IRFs is that in a nonlinear model, they will depend of the state of the model when the shock takes place and the expected value of the trejectory back to equilibrium is a nonlinear function of future shocks.
So Dynare computes an average Impulse Response Function in the following way. Assume, you want a 40 period IRF.
1) draw one history of shocks for 140 periods. Make two copies of it
2) In the second copy, add an extra 1 std err shock in period 101 on the exogenous variable that you want to compute IRF for.
3) Run a dynamic simulation with the two shock histories, starting at steady state
4) Compute the difference between the two trajectories of each endogenous variables
5) Repeat 1)-4) 50 times (Dynare default for option replic) and average the individual IRFs
Remarks:
a) the first 100 periods before the deterministic 1 std err shock do that the shock takes place at a random state for the economy
b) the next 40 periods shocks give one possible history after the deterministic shock
c) averaging over 50 different histories provide an average IRFs over state and future shocks (after the deterministic impulsion)
This explains why the IRFs may not statisfy the equation of the model.
Kind regards
Michel