| OPTIMIZE {bimets} | R Documentation |
Optimal control of a BIMETS model
Description
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:
1) check the correctness of input arguments;
2) perform a STOCHSIMULATE by uniformly perturbing the INSTRUMENT variables inside the user-boundaries provided in the OptimizeBounds function argument;
3) 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;
4) 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:
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;
ii) The INSTRUMENT variables 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);
iii) The following restrictions are applied to the INSTRUMENT: g + cn^2/2 < 27 \wedge g + cn > 17, given cn as the Consumption "booster" (i.e. the add-factor) and g 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 g<20 due to the square root operation in the objective function, given instrument g \in (15,25).
Usage
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,
...)
Arguments
model |
see |
simAlgo |
see |
TSRANGE |
see |
simType |
see |
simConvergence |
see |
simIterLimit |
see |
ZeroErrorAC |
see |
BackFill |
see |
Exogenize |
see |
ConstantAdjustment |
see |
verbose |
see |
verboseSincePeriod |
see |
verboseVars |
see |
StochReplica |
see |
StochSeed |
see |
OptimizeBounds |
the named |
OptimizeRestrictions |
the named |
OptimizeFunctions |
the named |
quietly |
see |
RESCHECKeqList |
see |
JACOBIAN_SHOCK |
see |
JacobianDrop |
see |
forceForwardLooking |
see |
avoidCompliance |
see |
... |
see |
Value
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.
See Also
MDL
LOAD_MODEL
ESTIMATE
STOCHSIMULATE
MULTMATRIX
RENORM
TIMESERIES
BIMETS indexing
BIMETS configuration
Examples
#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