How theoretical moments are calculated?

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.ghuoo_.dr.ghu’

Is this right?

Almost. Please see [What are "theoretical moments" referring to?) and get back if you have additional questions.

The formula above is wrong.
In my understanding, it has to be changed as follows.

Let
X_t = A vector of all model variables
and
S_t_1 = A vector of predetermined state variables
and
U_t = A vector of exogenous shocks

Also,
Let oo_.dr.ghx = A and oo_.dr.ghu = C

dynare results is nothing but

X_t = A * S_t_1 + C * U_t

Now, unconditional variance of X_t is

A * Var[S_t_1] * A’ + C * Var[U_t] * C’

Since Var[U_t] is the covariance matrix (Sigma) of which each diagonal elements are those I assigned as magnitudes of shocks in the model,
above can be written as

A * Var[S_t_1] * A’ + C * Sigma * C’

However, I do not know how to handle Var[S_t_1].
If I replace Var[S_t_1] with Sigma as well,
computed theoretical variances are slightly more than those displayed in dynare results.

Is there something wrong above?

Dynare first computes the state transition matrices

[A,B] = kalman_transition_matrix(dr,ipred,1:nx,M_.exo_nbr);
It then solves the associated Lyapunov equation to get the variance of the state variances.

[vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,],],options_.debug);
Then defining

aa = ghx(iky,:); bb = ghu(iky,:);
you can compute the variance as

v(stationary_vars,stationary_vars) = aa*vx*aa'+ bb*M_.Sigma_e*bb';
At second order, things are a bit more complicated because of the constant term entering. The computations you are looking for are conducted in th_autocovariances.m