My code is as follows.

```
var y R pi z deltaq y_bar y_star deltae pi_star pih mc RDM PIDM ;
varexo epsR epsq epsy_star epspi_star epsz ;
parameters tau alfa rhoz kappa rhoR psi1 psi2 rhoq rhoy_star rhopi_star rss;
psi1 = 1.746 ;
psi2 = 0.536 ;
rhoR = 0.833 ;
alfa = 0.174 ;
rss = 2.545 ;
kappa = 1.055 ;
tau = 0.454 ;
rhoq = 0.274 ;
rhoz = 0.105 ;
rhopi_star = 0.437 ;
rhoy_star = 0.915 ;
sigma_q = 2.092 ;
sigma_z = 1.211 ;
sigma_ps = 2.779 ;
sigma_ys = 0.816 ;
sigma_r = 0.148 ;
model(linear);
// 1. IS-equation
y-y_bar = y(+1)-y_bar(+1)-(tau+alfa*(2-alfa)*(1-tau))*(R-pih(+1)-z(+1)) ;
// 2. Natural Level of Output
y_bar=-(alfa*(2-alfa)*(1-tau))/tau*y_star;
// 3. Philips Curve
pih = exp(-rss/400)*pih(+1) + kappa*mc ;
// 4. Marginal Cost
mc = (1/(tau+alfa*(2-alfa)*(1-tau)))*(y-y_bar) ;
// 5. Relationship between CPI and PPI
pi = pih-alfa*deltaq ;
// 6. Exchange rate
pi=deltae+(1-alfa)*deltaq+pi_star;
// 7. Monetary Policy Rule 1(a) : Baseline
R=rhoR*R(-1)+(1-rhoR)*(psi1*pi+psi2*y)+epsR;
// 8. Terms of Trade
deltaq=rhoq*deltaq(-1)+epsq;
// 9. Technology Difference
z=rhoz*z(-1)+epsz;
// 10. Foreign Output
y_star=rhoy_star*y_star(-1)+epsy_star;
// 11. Foreign Inflation
pi_star=rhopi_star*pi_star(-1)+epspi_star;
PIDM = 4*pi ;
RDM = 4*R ;
end;
steady;
check;
shocks ;
var epsq = sigma_q^2 ;
end ;
stoch_simul(irf=20,nograph) ;
```

Running the mod file above, I can see theoretical moments.

However, I do not understand how these results are computed.

In my understanding, theoretical variance is computed by

oo_.dr.ghu*[Covariance matrix of exogenous shock]*oo_.dr.ghu + oo_.dr.ghu*oo_.dr.ghu’

Is this right?