Simultaneous Equations Models in R I use the systemfit package to evaluate Simultaneous Equations Models (for example 2SLS). But I do not understand how to predict endogenous variables using the new values of exogenous variables using this package.

Therefore, I would like to know how to make prediction (including prediction intervals) using models of simultaneous equations in the "systemfit" package or in some other way (maybe someone knows other packages for this).

I am also interested in how I can calculate the elastic coefficients of models.

I will be grateful to the literature on the topic of modeling simultaneous equations in R.
If you search for “simultaneous equations” inside the Econometrics Task View web page you can find two results: the systemfit and the bimets packages.
As you stated, the systemfit package is a powerful tool for econometric estimation of simultaneous systems of linear and nonlinear equations, but it only provides fitting procedures, thus it cannot be used in your example in order to work out a forecast.
On the other hand, the bimets package implements, among others, simulation and forecasting procedures, and multiplier analysis; therefore you can easily calculate elasticities.
Take a look also at the Rstudio blog post about SEM in R:
https://rviews.rstudio.com/2021/01/22/sem-time-series-modeling/
An example of a simple Klein model in R follows:
#load library
library(bimets)

#define the Klein model    
kleinModelDef &lt;- &quot;
MODEL

COMMENT&gt; Modified Klein Model 1 of the U.S. Economy with PDL, 
COMMENT&gt; autocorrelation on errors, restrictions and conditional equation evaluations

COMMENT&gt; Consumption with autocorrelation on errors
BEHAVIORAL&gt; cn
TSRANGE 1923 1 1940 1
EQ&gt; cn =  a1 + a2*p + a3*TSLAG(p,1) + a4*(wp+wg) 
COEFF&gt; a1 a2 a3 a4
ERROR&gt; AUTO(2)

COMMENT&gt; Investment with restrictions
BEHAVIORAL&gt; i
TSRANGE 1923 1 1940 1
EQ&gt; i = b1 + b2*p + b3*TSLAG(p,1) + b4*TSLAG(k,1)
COEFF&gt; b1 b2 b3 b4
RESTRICT&gt; b2 + b3 = 1

COMMENT&gt; Demand for Labor with PDL
BEHAVIORAL&gt; wp 
TSRANGE 1923 1 1940 1
EQ&gt; wp = c1 + c2*(y+t-wg) + c3*TSLAG(y+t-wg,1) + c4*time
COEFF&gt; c1 c2 c3 c4
PDL&gt; c3 1 2

COMMENT&gt; Gross National Product
IDENTITY&gt; y
EQ&gt; y = cn + i + g - t

COMMENT&gt; Profits
IDENTITY&gt; p
EQ&gt; p = y - (wp+wg)

COMMENT&gt; Capital Stock with IF switches
IDENTITY&gt; k
EQ&gt; k = TSLAG(k,1) + i
IF&gt; i &gt; 0
IDENTITY&gt; k
EQ&gt; k = TSLAG(k,1) 
IF&gt; i &lt;= 0

END
&quot;

#load the model
kleinModel &lt;- LOAD_MODEL(modelText = kleinModelDef)

#define data
kleinModelData &lt;- list(  
  cn  =TIMESERIES(39.8,41.9,45,49.2,50.6,52.6,55.1,56.2,57.3,57.8,
                  55,50.9,45.6,46.5,48.7,51.3,57.7,58.7,57.5,61.6,65,69.7,  
                  START=c(1920,1),FREQ=1),
  g   =TIMESERIES(4.6,6.6,6.1,5.7,6.6,6.5,6.6,7.6,7.9,8.1,9.4,10.7,
                  10.2,9.3,10,10.5,10.3,11,13,14.4,15.4,22.3,   
                  START=c(1920,1),FREQ=1),
  i   =TIMESERIES(2.7,-.2,1.9,5.2,3,5.1,5.6,4.2,3,5.1,1,-3.4,-6.2,
                  -5.1,-3,-1.3,2.1,2,-1.9,1.3,3.3,4.9,  
                  START=c(1920,1),FREQ=1),
  k   =TIMESERIES(182.8,182.6,184.5,189.7,192.7,197.8,203.4,207.6,
                  210.6,215.7,216.7,213.3,207.1,202,199,197.7,199.8,
                  201.8,199.9,201.2,204.5,209.4,    
                  START=c(1920,1),FREQ=1),
  p   =TIMESERIES(12.7,12.4,16.9,18.4,19.4,20.1,19.6,19.8,21.1,21.7,
                  15.6,11.4,7,11.2,12.3,14,17.6,17.3,15.3,19,21.1,23.5, 
                  START=c(1920,1),FREQ=1),
  wp  =TIMESERIES(28.8,25.5,29.3,34.1,33.9,35.4,37.4,37.9,39.2,41.3,
                  37.9,34.5,29,28.5,30.6,33.2,36.8,41,38.2,41.6,45,53.3,    
                  START=c(1920,1),FREQ=1),
  y   =TIMESERIES(43.7,40.6,49.1,55.4,56.4,58.7,60.3,61.3,64,67,57.7,
                  50.7,41.3,45.3,48.9,53.3,61.8,65,61.2,68.4,74.1,85.3, 
                  START=c(1920,1),FREQ=1),
  t   =TIMESERIES(3.4,7.7,3.9,4.7,3.8,5.5,7,6.7,4.2,4,7.7,7.5,8.3,5.4,
                  6.8,7.2,8.3,6.7,7.4,8.9,9.6,11.6, 
                  START=c(1920,1),FREQ=1),
  time=TIMESERIES(NA,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,
                  1,2,3,4,5,6,7,8,9,10, 
                  START=c(1920,1),FREQ=1),
  wg  =TIMESERIES(2.2,2.7,2.9,2.9,3.1,3.2,3.3,3.6,3.7,4,4.2,4.8,
                  5.3,5.6,6,6.1,7.4,6.7,7.7,7.8,8,8.5,  
                  START=c(1920,1),FREQ=1)
);

#load time series into the model object
kleinModel &lt;- LOAD_MODEL_DATA(kleinModel,kleinModelData)
## Load model data &quot;kleinModelData&quot; into model &quot;kleinModelDef&quot;...
## ...LOAD MODEL DATA OK

#estimate the model
kleinModel &lt;- ESTIMATE(kleinModel)

#In order to forecast the model up to 1944, 
#we need to extend exogenous variables up to 1944
kleinModel$modelData &lt;- within(kleinModel$modelData,{
    wg    = TSEXTEND(wg,  UPTO=c(1944,1),EXTMODE='CONSTANT')
    t     = TSEXTEND(t,   UPTO=c(1944,1),EXTMODE='LINEAR')
    g     = TSEXTEND(g,   UPTO=c(1944,1),EXTMODE='CONSTANT')
    k     = TSEXTEND(k,   UPTO=c(1944,1),EXTMODE='LINEAR')
    time  = TSEXTEND(time,UPTO=c(1944,1),EXTMODE='LINEAR')
  })

#forecast model
kleinModel &lt;- SIMULATE(kleinModel
                      ,simType='FORECAST'
                      ,TSRANGE=c(1941,1,1944,1)
                      ,simConvergence=0.00001
                      ,simIterLimit=100
                      ,quietly=TRUE
  ) 

#get forecasted GNP
TABIT(kleinModel$simulation$y)
## 
##       DATE, PER, kleinModel$simulation$y
## 
##       1941, 1  ,  125.3      
##       1942, 1  ,  172.5      
##       1943, 1  ,  185.6      
##       1944, 1  ,  141.1