Here’s the new plan:
For each parameter, solve for the state space form (z_t,g_t,R_t,y_t,\pi_t) '= A\cdot (z_{t-1},g_{t-1},R_{t-1},y_{t-1},\pi_{t-1}) '+B\cdot (e_{t}^z,e_{t}^g,e_{t}^R)'.
E_{t=0}(y_{t+s}^2+\pi_{t+s}^2)= \left[A(4,:)\cdot A(1:3,:)^{s-1}\cdot (z_{t=0},g_{t=0},R_{t=0})'\right]^2+ \left[A(5,:)\cdot A(1:3,:)^{s-1}\cdot (z_{t=0},g_{t=0},R_{t=0})'\right]^2
+ \text{ some term that consists of } \sigma_z, \sigma_R, \sigma_g
No need to simulate anything, just solve and compute. Sounds reasonable?