# SEM Stochastic Simulation and Optimal Control

**R Views**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

*Andrea Luciani is a Technical Advisor for the Directorate General for Economics, Statistics and Research at the Bank of Italy, and co-author of the bimets package.*

#### Introduction

In this post, I show how to analyze the forecast error for a Structural Equation Model (SEM) by means of a stochastic simulation, and how to perform optimal control. In my previous post, I presented an approach to estimating and simulating SEMs in R which focused on R tools that allow users to forecast advanced simultaneous equations models having linear restrictions on coefficients, error autocorrelation on residuals, and conditional equation evaluation. I described how simulating these econometric models often requires using iterative algorithms in ways that are well beyond what the `ts()`

, `lm()`

and `predict()`

functions were designed to do and worked through a forecasting exercise of a Klein model by using the deterministic simulation procedure available into the bimets package.

#### Stochastic Simulation

Deterministic algorithms, however, have significant limitations when applied to SEM forecasts which are subject to several sources of error such as: a random disturbance term of each stochastic equation, errors in estimated coefficients, and errors in forecasts of exogenous variables.

The forecast error depending on the structural disturbances can be analyzed by using a stochastic simulation in which the structural disturbances are given values with specified stochastic properties. The error terms of the estimated equations are appropriately perturbed. Technical identities and exogenous variables can also be perturbed by disturbances with specified stochastic properties. Programmatically, a straightforward approach to solve the problem is a Monte-Carlo method: the model is solved for each data set with different values for the disturbances. Finally, forecast statistics (e.g. mean, standard deviation, etc.) can be computed for each simulated endogenous variable.

Using the Klein model defined in my previous post, we can attempt to complete a US *Gross National Product* forecasting exercise, given a user-specified uncertainty on *Consumption* and *Government Expenditure* time series.

The Klein model will be perturbed by applying a normal disturbance to the endogenous *Consumption* behavioral `cn`

in year 1942, and a uniform disturbance to the exogenous *Government Expenditure* time series `g`

along all the forecast time range. The normal disturbance applied to the `cn`

stochastic equation has a zero mean and a standard deviation equal to its regression standard error, thus roughly replicating the regression error during the current perturbation (not accounting for inter-equations cross-covariance).

#we want to perform a stochastic forecast of the GNP up to 1944 #we will add normal disturbances to endogenous Consumption 'cn' #in 1942 by using its regression standard error #we will add uniform disturbance to exogenous Government Expenditure 'g' #in whole TSRANGE myStochStructure <- list( cn=list( TSRANGE=c(1942,1,1942,1), TYPE='NORM', PARS=c(0,kleinModel$behaviorals$cn$statistics$StandardErrorRegression) ), g=list( TSRANGE=TRUE, TYPE='UNIF', PARS=c(-1,1) ) ) #model stochastic forecast kleinModel <- STOCHSIMULATE(kleinModel ,simType='FORECAST' ,TSRANGE=c(1941,1,1944,1) ,StochStructure=myStochStructure ,StochSeed=123 ,quietly=TRUE ) #print mean and standard deviation for the forecasted GNP with(kleinModel$stochastic_simulation,TABIT(y$mean, y$sd)) ## ## Date, Prd., y$mean , y$sd ## ## 1941, 1 , 125.5 , 4.251 ## 1942, 1 , 173.3 , 9.263 ## 1943, 1 , 186 , 11.88 ## 1944, 1 , 141.1 , 11.7

#### Optimal Control

Another limitation of the deterministic simulation procedure relies on the fact that analysts often do not want to perform a mere forecasting exercise. Usually, a performance measure (i.e. *objective-function*) is assigned to an econometric model, depending on the value of forecasted endogenous variables; thus, analysts try to enhance this measure by fine-tuning exogenous variables (i.e. the *instruments*) of the model (e.g. policy analysis).

Generally speaking, this is an *optimal control* exercise: the optimization consists of maximizing an objective-function, depending on (forecasted) endogenous variables, given a set of user-constraints plus the constraints imposed by the econometric model equations.

In the following exercise, we will use the bimets `OPTIMIZE()`

procedure which allows R users to perform optimal control exercises on simultaneous equations models. The exercise will be completed by using the following code on the previously defined Klein model, and we will assume that:

The objective-function definition is: \(f(y, cn, g) = (y-110)+(cn-90)*|cn-90|-\sqrt{g-20}\) given \(y\) as the simulated

*Gross National Product*, \(cn\) as the simulated*Consumption*and \(g\) as the exogenous*Government Expenditure*. The basic idea is to maximize*Consumption*, and secondarily the*Gross National Product*, while reducing the*Government Expenditure*;The instrument variables (i.e.

`INSTRUMENT`

in following code snippet) are the \(cn\)*Consumption*“booster” (i.e. the*add-factor*, not to be confused with the simulated*Consumption*in the objective-function) and the \(g\)*Government Expenditure*, defined over the following domains: \(cn \in (-5,5)\), \(g \in (15,25)\);The following restrictions are applied to the

`INSTRUMENT`

: \(g + cn^2/2 < 27 \wedge g + cn > 17\), given again \(cn\) as the*Consumption*“booster” (i.e. the*add-factor*) and \(g\) as the*Government Expenditure*;

The following code performs the optimization.

#we want to maximize the non-linear objective function: #f()=(y-110)+(cn-90)*ABS(cn-90)-(g-20)^0.5 #in 1942 by using INSTRUMENT cn in range (-5,5) #(cn is endogenous so we use the add-factor) #and g in range (15,25) #we will also impose the following non-linear restriction: #g+(cn^2)/2<27 & g+cn>17 #define INSTRUMENT and boundaries myOptimizeBounds <- list( cn = list( TSRANGE = TRUE ,BOUNDS = c(-5,5)), g = list( TSRANGE = TRUE ,BOUNDS = c(15,25)) ) #define restrictions myOptimizeRestrictions <- list( myRes1=list( TSRANGE = TRUE ,INEQUALITY = 'g+(cn^2)/2<27 & g+cn>17') ) #define objective function myOptimizeFunctions <- list( myFun1 = list( TSRANGE = TRUE ,FUNCTION = '(y-110)+(cn-90)*ABS(cn-90)-(g-20)^0.5') ) #Monte-Carlo optimization by using 10000 stochastic realizations #and 1E-4 convergence criterion kleinModel <- OPTIMIZE(kleinModel ,simType = 'FORECAST' ,TSRANGE=c(1942,1,1942,1) ,simConvergence= 1E-4 ,simIterLimit = 1000 ,StochReplica = 10000 ,StochSeed = 123 ,OptimizeBounds = myOptimizeBounds ,OptimizeRestrictions = myOptimizeRestrictions ,OptimizeFunctions = myOptimizeFunctions ) ## OPTIMIZE(): optimization boundaries for the add-factor of endogenous variable "cn" are (-5,5) from year-period 1942-1 to 1942-1. ## OPTIMIZE(): optimization boundaries for the exogenous variable "g" are (15,25) from year-period 1942-1 to 1942-1. ## OPTIMIZE(): optimization restriction "myRes1" is active from year-period 1942-1 to 1942-1. ## OPTIMIZE(): optimization objective function "myFun1" is active from year-period 1942-1 to 1942-1. ## ## Optimize: 100.00 % ## OPTIMIZE(): 2916 out of 10000 objective function realizations (29%) are finite and verify the provided restrictions. ## ...OPTIMIZE OK #print local maximum kleinModel$optimize$optFunMax ## [1] 210.6 #print INSTRUMENT that allow local maximum to be achieved kleinModel$optimize$INSTRUMENT ## $cn ## Time Series: ## Start = 1942 ## End = 1942 ## Frequency = 1 ## [1] 2.032 ## ## $g ## Time Series: ## Start = 1942 ## End = 1942 ## Frequency = 1 ## [1] 24.9

**leave a comment**for the author, please follow the link and comment on their blog:

**R Views**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.