Optimal control of a BIMETS model
The OPTIMIZE procedure provides a convenient method for performing optimal control exercises; the procedure maximizes an arbitrary objective-function under the constraints imposed by the econometric model and by user-specified constraints.
An approach to policy evaluation is via a so-called "social welfare function". This approach relaxes the assumptions of the instruments-targets framework, i.e. the RENORM procedure. Rather than assuming specific desired targets for some endogenous variables, it assumes the existence of a social welfare function determining a scalar measure of performance based on both endogenous and policy (exogenous) variables.
The social welfare function can incorporate information about tradeoffs in objectives that are not allowed by the RENORM instruments-targets approach.
BIMETS supplies the OPTIMIZE procedure in order to perform optimal control exercises on econometric models.
The optimization consists of maximizing a social welfare function, i.e. the objective-function, depending on exogenous and (simulated) endogenous variables, subject to user constraints plus the constraints imposed by the econometric model equations. Users are allowed to define constraints and objective-functions of any degree, and are allowed to provide different constraints and objective-functions in different optimization time periods.
The core of the OPTIMIZE procedure is based on a Monte Carlo method that takes advantage of the STOCHSIMULATE procedure. Policy variables, i.e. INSTRUMENT, are uniformly perturbed in the range defined by the user-provided boundaries, then the INSTRUMENT values that i) verify the user-provided constraints and ii) maximize the objective-functions are selected and stored into the optimize element of the output BIMETS model.
The following steps can describe the procedure implemented in OPTIMIZE:
check the correctness of input arguments;
perform a STOCHSIMULATE by uniformly perturbing the INSTRUMENT variables inside the user-boundaries provided in the OptimizeBounds function argument;
during the STOCHSIMULATE, for each period in the optimization TSRANGE: i) discard the stochastic realizations that do not verify the restrictions provided in the OptimizeRestrictions argument; ii) for all the remaining realizations, compute the current value of the objective-functions time series, as defined in the OptimizeFunctions argument, by using the exogenous and (simulated) endogenous stochastic time series;
once the STOCHSIMULATE completes, select the stochastic realization that presents the higher value in the sum of the corresponding objective-function time series values, and return, among other data, the related optimal INSTRUMENT time series.
In the following figure, the scatter plot is populated with 2916 objective function stochastic realizations, computed by using the example code at the end of this section; the 210.58 local maximum is highlighted
(i.e. advancedKleinModel$optimize$optFunMax in first example).
In this example:
i) The objective function definition is:
given as the simulated Gross National Product, as the simulated Consumption and as the exogenous Government Expenditure: the basic idea is to maximize Consumption, and secondarily the Gross National Product, while reducing the Government Expenditure;
ii) The INSTRUMENT variables are the Consumption "booster" (i.e. the add-factor, not to be confused with the simulated Consumption in the objective function) and the Government Expenditure, defined over the following domains: , ;
iii) The following restrictions are applied to the INSTRUMENT: , given as the Consumption "booster" (i.e. the add-factor) and as the Government Expenditure;
The figure clearly shows that non-linear restrictions have been applied, and that non-computable objective functions have been discarded, e.g. the stochastic realizations having due to the square root operation in the objective function, given instrument .
OPTIMIZE( model=NULL, simAlgo='GAUSS-SEIDEL', TSRANGE=NULL, simType='DYNAMIC', simConvergence=0.01, simIterLimit=100, ZeroErrorAC=FALSE, BackFill=0, Exogenize=NULL, ConstantAdjustment=NULL, verbose=FALSE, verboseSincePeriod=0, verboseVars=NULL, StochReplica=100, StochSeed=NULL, OptimizeBounds=NULL, OptimizeRestrictions=NULL, OptimizeFunctions=NULL, quietly=FALSE, RESCHECKeqList=NULL, JACOBIAN_SHOCK=1e-4, JacobianDrop=NULL, forceForwardLooking=FALSE, avoidCompliance=FALSE, ...)
model: see SIMULATE
simAlgo: see SIMULATE
TSRANGE: see SIMULATE
simType: see SIMULATE
simConvergence: see SIMULATE
simIterLimit: see SIMULATE
ZeroErrorAC: see SIMULATE
BackFill: see SIMULATE
Exogenize: see SIMULATE
ConstantAdjustment: see SIMULATE
verbose: see SIMULATE
verboseSincePeriod: see SIMULATE
verboseVars: see SIMULATE
StochReplica: see STOCHSIMULATE
StochSeed: see STOCHSIMULATE
OptimizeBounds: the named list() that defines the search boundaries applied to INSTRUMENT exogenous variables. Each list element must have a name equal to an endogenous or an exogenous model variable.
The list names define the INSTRUMENT.
If a list element name is equal to an exogenous variable, then the boundaries will be applied directly to the related exogenous stochastic time series values. If a list element name is equal to an endogenous variable, then the boundaries will be applied to the stochastic constant adjustment (see STOCHSIMULATE) of the related endogenous variable.
Each list element must be a named list built with the following two named variables:
TSRANGE: the time range wherein the search boundaries are active. The TSRANGE must be a 4 numerical array,i.e. TSRANGE=c(start_year, start_period, end_year, end_period) or TSRANGE=TRUE in order to apply the provided boundaries to the whole OPTIMIZE TSRANGE.
BOUNDS: the boundaries that are applied to the related instrument. These parameters must contain the lower and upper bound of the uniform distribution wherein the search for the objective-functions maximum is performed,i.e. BOUNDS=c(lower_bound,upper_bound).
See example in order to learn how to build a compliant boundaries structure.
OptimizeRestrictions: the named list() that defines the restrictions applied to INSTRUMENT exogenous variables. This list can be NULL.
Each list element must be a named list built with the following two named variables:
TSRANGE: the time range wherein the restriction is active. The TSRANGE must be a 4 numerical array,i.e. TSRANGE=c(start_year, start_period, end_year, end_period) or TSRANGE=TRUE in order to apply the provided restriction to the whole OPTIMIZE TSRANGE.
INEQUALITY: the inequality expression, i.e. a character variable, that defines the restriction. The INEQUALITY expression can contain exogenous and endogenous variable names, the standard arithmetic and logical operators, parentheses and the MDL functions described in the EQ section of the MDL help page. If in the INEQUALITY expression a variable name refers to an exogenous variable, then that variable will be evaluated by using the related exogenous time series stochastic values. If in the INEQUALITY expression a variable name refers to an endogenous variable, then that variable will be evaluated by using to the stochastic constant adjustment (see argument StochStructure of the STOCHSIMULATE help page) of the related endogenous variable.Two different OptimizeRestrictions list element can not have overlapping TSRANGE. See example in order to learn how to build a compliant restrictions structure.
OptimizeFunctions: the named list() that defines the objective functions to be maximized.
Each list element must be a named list built with the following two named variables:
TSRANGE: the time range wherein the objective function is evaluated. The TSRANGE must be a 4 numerical array,i.e. TSRANGE=c(start_year, start_period, end_year, end_period) or TSRANGE=TRUE in order to evaluate the objective function in each period of the OPTIMIZE TSRANGE.
FUNCTION: the expression, i.e. a character variable, that defines the objective function. The FUNCTION expression can contain exogenous and endogenous variable names, the standard arithmetic and logical operators, parentheses and the MDL functions described in the EQ section of the MDL help page. If in the FUNCTION expression a variable name refers to an exogenous variable, then that variable will be evaluated by using the related exogenous time series stochastic values. If in the FUNCTION expression a variable name refers to an endogenous variable, then that variable will be evaluated by using the stochastic simulated time series (see STOCHSIMULATE) of the related endogenous variable.Two different OptimizeFunctions list element can not have overlapping TSRANGE. See example in order to learn how to build a compliant objective functions structure.
quietly: see SIMULATE
RESCHECKeqList: see SIMULATE
JACOBIAN_SHOCK: see SIMULATE
JacobianDrop: see SIMULATE
forceForwardLooking: see SIMULATE
avoidCompliance: see SIMULATE
...: see SIMULATE
This function will add, into the output BIMETS model object, three new named elements, respectively optimize, simulation_MM and INSTRUMENT_MM.
The optimize element is a named list() that contains the following elements:
INSTRUMENT: a named list that contains the time series of the instrument exogenous variables that verify the OptimizeRestrictions and that allow the objective OptimizeFunctions to be maximized. This element is populated only if a finite solution exists. List names are equal to the names of the related exogenous variables. Users can also declare an endogenous variable as INSTRUMENT variable, by using the OptimizeBounds argument: in this case the constant adjustment (see STOCHSIMULATE) related to the provided endogenous variable will be used as instrument exogenous variable, and this output INSTRUMENT list will contains the constant adjustment time series that allow the objective OptimizeFunction to be maximized (see example);
optFunMax: the scalar value (local maximum) obtained by evaluating the OptimizeFunctions while the model is fed by the optimized INSTRUMENT time series. This element is populated only if a finite solution exists;
optFunTS: the time series obtained by evaluating the OptimizeFunctions during each period in the OPTIMIZE TSRANGE while the model is fed by the optimized INSTRUMENT time series. Thus, optFunMax==sum(optFunTS). This element is populated only if a finite solution exists;
optFunAve: the scalar value that is the mean of all the stochastic OptimizeFunctions realizations, filtered by the restrictions imposed by the OptimizeRestrictions argument. This element is populated only if a finite solution exists;
optFunSd: the scalar value that is the standard deviation of all the stochastic OptimizeFunctions realizations, filtered by the restrictions imposed by the OptimizeRestrictions argument. This element is populated only if a finite solution exists;
realizationsToKeep: a 1 x StochReplica boolean row array. If the i-th element is TRUE than the related objective function realization is computable and verifies the restrictions imposed by the OptimizeRestricions argument. It can be useful along with optFunResults and INSTRUMENT_MM in order to verify and to refine results;
optFunResults: the numerical array containing the evaluated OptimizeFunctions for all the (unfiltered) realizations;
modelData: the whole model input dataset wherein the INSTRUMENT exogenous variables have been modified accordingly to the OPTIMIZE results. This data can be useful in order to verify or to refine results (see example);
ConstantAdjustment: a modified constant adjustment input list wherein the constant adjustment time series related to a INSTRUMENT endogenous variables have been modified accordingly to the OPTIMIZE results. This data can be useful in order to verify or to refine results (see example);
The arguments passed to the function call during the latest OPTIMIZE run will be inserted into the '__OPT_PARAMETERS__' element of the model optimize list; this data can be helpful in order to replicate the optimization results.
The simulation_MM element is a named list(), having the endogenous variables as names. Each element will contain an R x C matrix, given R the number of observations in the optimization TSRANGE and C=1+StochReplica. The first column of each matrix contains the related endogenous variable's unperturbed simulated values; the remaining columns will contain all the StochReplica stochastic realizations for the related endogenous variable.
The INSTRUMENT_MM element is a named list(), having INSTRUMENT variables as names. Each element will contain an R x C matrix, given R the number of observations in the optimization TSRANGE and C=1+StochReplica. The first column of each matrix contains the related INSTRUMENT variable's unperturbed values; the remaining columns will contain all the StochReplica stochastic realizations for the related INSTRUMENT variable.
MDL
LOAD_MODEL
ESTIMATE
STOCHSIMULATE
MULTMATRIX
RENORM
TIMESERIES
BIMETS indexing
BIMETS configuration
#define the advanced Klein model advancedKleinModelDef <- " MODEL COMMENT> Modified Klein Model 1 of the U.S. Economy with PDL, COMMENT> autocorrelation on errors, restrictions and conditional equation evaluations COMMENT> Consumption with autocorrelation on errors BEHAVIORAL> cn TSRANGE 1923 1 1940 1 EQ> cn = a1 + a2*p + a3*TSLAG(p,1) + a4*(w1+w2) COEFF> a1 a2 a3 a4 ERROR> AUTO(2) COMMENT> Investment with restrictions BEHAVIORAL> i TSRANGE 1923 1 1940 1 EQ> i = b1 + b2*p + b3*TSLAG(p,1) + b4*TSLAG(k,1) COEFF> b1 b2 b3 b4 RESTRICT> b2 + b3 = 1 COMMENT> Demand for Labor with PDL BEHAVIORAL> w1 TSRANGE 1923 1 1940 1 EQ> w1 = c1 + c2*(y+t-w2) + c3*TSLAG(y+t-w2,1) + c4*time COEFF> c1 c2 c3 c4 PDL> c3 1 2 COMMENT> Gross National Product IDENTITY> y EQ> y = cn + i + g - t COMMENT> Profits IDENTITY> p EQ> p = y - (w1+w2) COMMENT> Capital Stock with IF switches IDENTITY> k EQ> k = TSLAG(k,1) + i IF> i > 0 IDENTITY> k EQ> k = TSLAG(k,1) IF> i <= 0 END " #load the model advancedKleinModel <- LOAD_MODEL(modelText = advancedKleinModelDef) #define data kleinModelData <- 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), w1 =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), w2 =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 advancedKleinModel <- LOAD_MODEL_DATA(advancedKleinModel,kleinModelData) #estimate the model advancedKleinModel <- ESTIMATE(advancedKleinModel, quietly=TRUE) #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 #we need to extend exogenous variables up to 1942 advancedKleinModel$modelData <- within(advancedKleinModel$modelData,{ w2 = TSEXTEND(w2, UPTO=c(1942,1),EXTMODE='CONSTANT') t = TSEXTEND(t, UPTO=c(1942,1),EXTMODE='LINEAR') g = TSEXTEND(g, UPTO=c(1942,1),EXTMODE='CONSTANT') k = TSEXTEND(k, UPTO=c(1942,1),EXTMODE='LINEAR') time = TSEXTEND(time,UPTO=c(1942,1),EXTMODE='LINEAR') }) #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 advancedKleinModel <- OPTIMIZE(advancedKleinModel ,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 advancedKleinModel$optimize$optFunMax #[1] 210.5755 #print INSTRUMENT that allow local maximum to be achieved advancedKleinModel$optimize$INSTRUMENT #$cn #Time Series: #Start = 1942 #End = 1942 #Frequency = 1 #[1] 2.032203 # #$g #Time Series: #Start = 1942 #End = 1942 #Frequency = 1 #[1] 24.89773 #LET'S VERIFY RESULTS #copy into modelData the computed INSTRUMENT #that allow to maximize the objective function advancedKleinModel$modelData <- advancedKleinModel$optimize$modelData #simulate the model by using the new INSTRUMENT #note: we used cn add-factor as OPTIMIZE instrument, so we need #to pass the computed cn add-factor to the SIMULATE call newConstantAdjustment <- advancedKleinModel$optimize$ConstantAdjustment advancedKleinModel <- SIMULATE(advancedKleinModel ,simType = 'FORECAST' ,TSRANGE = c(1942,1,1942,1) ,simConvergence = 1E-5 ,simIterLimit = 1000 ,ConstantAdjustment = newConstantAdjustment ) #calculate objective function by using the SIMULATE output time series #(y-110)+(cn-90)*ABS(cn-90)-(g-20)^0.5 y <- advancedKleinModel$simulation$y cn <- advancedKleinModel$simulation$cn g <- advancedKleinModel$modelData$g optFunTest <- (y-110)+(cn-90)*abs(cn-90)-(g-20)^0.5 #verify computed max is equal to optimization max #(in the following command TSPROJECT could be omitted because #myFun1$TSRANGE = TRUE) abs(sum(TSPROJECT(optFunTest ,TSRANGE=c(1942,1,1942,1) ,ARRAY = TRUE) ) - advancedKleinModel$optimize$optFunMax) < 1E-4 #[1] TRUE #we can also check that the SIMULATE time series #are equal to the OPTIMIZE realizations that allow to maximize #the objective function #get realization index that maximizes the objective function maximizingRealizationIdx <- with(advancedKleinModel$optimize, which.max(optFunResults[realizationsToKeep])) #get stochastic realizations unfiltered #(simulation_MM and INSTRUMENT_MM are populated during the OPTIMIZE call) y_opt <- advancedKleinModel$simulation_MM$y cn_opt <- advancedKleinModel$simulation_MM$cn g_opt <- advancedKleinModel$INSTRUMENT_MM$g #filter by restrictions and by finite solutions #(first column in all matrices is related to the un-perturbed model) y_opt <- y_opt[ ,c(FALSE,advancedKleinModel$optimize$realizationsToKeep),drop=FALSE] cn_opt <- cn_opt[,c(FALSE,advancedKleinModel$optimize$realizationsToKeep),drop=FALSE] g_opt <- g_opt[ ,c(FALSE,advancedKleinModel$optimize$realizationsToKeep),drop=FALSE] #get maximizing realizations y_opt <- y_opt[ ,maximizingRealizationIdx,drop=FALSE] cn_opt <- cn_opt[,maximizingRealizationIdx,drop=FALSE] g_opt <- g_opt[ ,maximizingRealizationIdx,drop=FALSE] #verify that these variables are equal to the SIMULATE time series max(abs(y-y_opt)) < 1E-4 #[1] TRUE max(abs(cn-cn_opt)) < 1E-4 #[1] TRUE max(abs(g[[1942,1]]-g_opt)) < 1E-4 #[1] TRUE ############################################################ #MULTI RESTRICTIONS, MULTI OBJECTIVE FUNCTIONS EXAMPLE #load the model (reset stuff) advancedKleinModel <- LOAD_MODEL(modelText = advancedKleinModelDef) #load time series into the model object advancedKleinModel <- LOAD_MODEL_DATA(advancedKleinModel,kleinModelData) #estimate the model advancedKleinModel <- ESTIMATE(advancedKleinModel, quietly=TRUE) #we want to maximize the non-linear objective function: #f1()=(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 want to maximize the non-linear objective function: #f2()=(y-120)+(cn-100)*ABS(cn-100)-(g-20)^0.5-(w2-8)^0.5 #in 1943 by using INSTRUMENT cn in range (-5,5), #g in range (15,25) #and w2 in range (7.5,12.5) #we will also impose the following non-linear restrictions: #in 1942: g+(cn^2)/2<27 & g+cn>17 #in 1943: (g^2)/10+(cn^2)/2+w2^2 < 200 #we need to extend exogenous variables up to 1943 advancedKleinModel$modelData <- within(advancedKleinModel$modelData,{ w2 = TSEXTEND(w2, UPTO=c(1943,1),EXTMODE='CONSTANT') t = TSEXTEND(t, UPTO=c(1943,1),EXTMODE='LINEAR') g = TSEXTEND(g, UPTO=c(1943,1),EXTMODE='CONSTANT') k = TSEXTEND(k, UPTO=c(1943,1),EXTMODE='LINEAR') time = TSEXTEND(time,UPTO=c(1943,1),EXTMODE='LINEAR') }) #define INSTRUMENT and boundaries myOptimizeBounds <- list( cn=list(TSRANGE=TRUE, BOUNDS=c(-5,5)), g=list(TSRANGE=TRUE, BOUNDS=c(15,25)), w2=list(TSRANGE=c(1943,1,1943,1), BOUNDS=c(7.5,12.5)) ) #define restrictions myOptimizeRestrictions <- list( myRes1=list( TSRANGE=c(1942,1,1942,1), INEQUALITY='g+(cn^2)/2 < 27 & g+cn > 17'), myRes2=list( TSRANGE=c(1943,1,1943,1), INEQUALITY='(g^2)/10+(cn^2)/2+w2^2 < 200') ) #define objective functions myOptimizeFunctions <- list( myFun1=list( TSRANGE=c(1942,1,1942,1), FUNCTION='(y-110)+(cn-90)*ABS(cn-90)-(g-20)^0.5'), myFun2=list( TSRANGE=c(1943,1,1943,1), FUNCTION='(y-120)+(cn-100)*ABS(cn-100)-(g-20)^0.5-(w2-8)^0.5') ) #Monte-Carlo optimization by using 1000 stochastic realizations #and 1E-4 convergence advancedKleinModel <- OPTIMIZE(advancedKleinModel ,simType = 'FORECAST' ,TSRANGE=c(1942,1,1943,1) ,simConvergence=1E-4 ,simIterLimit = 500 ,StochReplica = 1000 ,StochSeed = 123 ,OptimizeBounds = myOptimizeBounds ,OptimizeRestrictions = myOptimizeRestrictions ,OptimizeFunctions = myOptimizeFunctions) #print INSTRUMENT that allow local maximum to be achieved advancedKleinModel$optimize$INSTRUMENT #LET'S VERIFY RESULTS #copy into modelData the computed INSTRUMENT #that allow to maximize the objective function advancedKleinModel$modelData <- advancedKleinModel$optimize$modelData #simulate the model by using the new INSTRUMENT newConstantAdjustment <- advancedKleinModel$optimize$ConstantAdjustment advancedKleinModel <- SIMULATE(advancedKleinModel ,simType = 'FORECAST' ,TSRANGE = c(1942,1,1943,1) ,simConvergence = 1E-5 ,simIterLimit = 100 ,ConstantAdjustment = newConstantAdjustment ) #calculate objective functions by using the SIMULATE output time series y <- advancedKleinModel$simulation$y cn <- advancedKleinModel$simulation$cn g <- advancedKleinModel$modelData$g w2 <- advancedKleinModel$modelData$w2 optFunTest1 <- (y-110)+(cn-90)*abs(cn-90)-(g-20)^0.5 optFunTest2 <- (y-120)+(cn-100)*abs(cn-100)-(g-20)^0.5-(w2-8)^0.5 #verify computed max is equal to optimization max abs(sum(TSPROJECT(optFunTest1 ,TSRANGE=c(1942,1,1942,1) ,ARRAY = TRUE)+ TSPROJECT(optFunTest2 ,TSRANGE=c(1943,1,1943,1) ,ARRAY = TRUE) ) - advancedKleinModel$optimize$optFunMax) < 1E-2 #[1] TRUE