SEM Stochastic Simulation and Optimal Control

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:

1. 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;

2. 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)\);

3. 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

The scatter plot of objective function surface is populated with 2916 objective-function stochastic realizations; the 210.58 local maximum is stored in the kleinModel$optimize$optFunMax variable, while the instruments that allow local maximum to be achieved are stored in the kleinModel$optimize$INSTRUMENT variable.

Summary

Forecast errors in Structural Equation Models can be analyzed by using a stochastic simulation in which the structural disturbances are given values with specified stochastic properties. Stochastic modeling forecasts the probability of various outcomes under different conditions, using random variables. Often, a performance measure (i.e. objective-function) is assigned to an econometric model, depending on the value of forecasted endogenous variables. The Monte Carlo simulation is one example of a stochastic model; it can simulate how an econometric model may perform based on the probability distributions of individual variables of the model, allowing analysts to enhance the objective-function by fine-tuning exogenous variables (i.e. the instruments) of the model (e.g. policy analysis) in the so-called optimal control procedure. Examples presented in this post show how to perform stochastic simulation and optimal control in R.

Disclaimer: The views and opinions expressed in this page are those of the author and do not necessarily reflect the official policy or position of the Bank of Italy. Examples of analysis performed within these pages are only examples. They should not be utilized in real-world analytic products as they are based only on very limited and dated open source information. Assumptions made within the analysis are not reflective of the position of the Bank of Italy.