STOCHSIMULATE {bimets}

Stochastic simulation of a BIMETS model

Description

The STOCHSIMULATE operation performs a stochastic simulation. The simulation algorithms are the same as those used by the SIMULATE operation.

Forecasts produced by structural econometric models are subject to several sources of error, such as random disturbance term of each stochastic equation, errors in estimated coefficients, errors in forecasts of exogenous variables, errors in preliminary data and mis-specification of the model.

The forecast error depending on the structural disturbances can be analyzed by using the stochastic simulation procedure.

The deterministic simulation is the simultaneous solution of an econometric model obtained by applying, for each stochastic (behavioral) equation, the expected values of the structural disturbances, which are all zero by assumption. In the BIMETS STOCHSIMULATE stochastic simulation, the structural disturbances are given values that have specified stochastic properties. The error terms of the estimated behavioral equation of the model are appropriately perturbed. Identity equations and exogenous variables can be as well perturbed by disturbances that have specified stochastic properties. The model is then solved for each data set with different values of the disturbances. Finally, mean and standard deviation are computed for each simulated endogenous variable.

In terms of computational efficiency, the procedure takes advantage of the fact that multiple datasets are bound together in matrices. Therefore, to achieve a global convergence, the iterative simulation algorithm is executed once for all perturbed datasets. This solution can be viewed as a sort of a SIMD (i.e. Single Instruction Multiple Data) parallel simulation: the STOCHSIMULATE function transforms time series into matrices and consequently can easily bind multiple datasets by column. At the same time, a single run ensures a fast code execution. Finally, each column in the output matrices represents a stochastic realization.

By using the StochStructure argument of this function, users can define a stochastic structure for the disturbances. For each variable of the model, users can provide a distinct probability distribution for the disturbance, and a specific time range of application. Mean and standard deviation for each simulated endogenous time series will be stored in the stochastic_simulation element of the output model object; all the stochastic realizations will be stored in the simulation_MM element of the output model object as named matrices.

In the following figure, the advanced Klein model (see code example), has been perturbed during the forecast operation by applying a normal disturbance to the endogenous Consumption behavioral cn add-factor in year 1942, and a uniform disturbance to the exogenous Government Expenditure time series g along all the simulation TSRANGE. The normal disturbance applied to the cn behavioral has a zero mean and a standard deviation equal to the behavioral regression standard error,
i.e. advancedKleinModel$behaviorals$cn$statistics$StandardErrorRegression, thus roughly replicating the ESTIMATE regression error during the current perturbation (not accounting for inter-equations cross-covariance).

At the moment, all the disturbances are i.i.d. and are not transformed into a congruent autoregressive scheme in the case the related perturbed endogenous behavioral presents an autocorrelation for the errors in its MDL definition, e.g. ERROR> AUTO(n)

For more realistic scenarios, several advanced econometric exercises on the US Federal Reserve FRB/US econometric model (e.g., dynamic simulation in a monetary policy shock, rational expectations, endogenous targeting, stochastic simulation, etc.) are available in the "US Federal Reserve quarterly model (FRB/US) in R with bimets" vignette.

Usage

STOCHSIMULATE(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,
              StochStructure=NULL,
              StochReplica=100,
              StochSeed=NULL,
              quietly=FALSE,
              RESCHECKeqList=NULL,
              JACOBIAN_SHOCK=1e-4,
              JacobianDrop=NULL,
              forceForwardLooking=FALSE,
              avoidCompliance=FALSE,
              ...)		   

Arguments

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

StochStructure

the named list() that defines the disturbance structure applied to the model. Each list element must have a name equal to an endogenous or an exogenous model variable. List names define the INSTRUMENT.

Generic syntax:

StochStructure = 

list(

  var_name1 = list(
      TSRANGE=...,
      TYPE=...,
      PARS=...),
      
  var_name2 = list(...)
  
  ...
  
  var_nameN = list(...)
  
  )
  
  

ENDOGENOUS REFERENCE.

If the generic name var_nameI refers to an endogenous variable, the disturbance will be additive and applied to the constant adjustment time series (see SIMULATE) of the referred endogenous variable.

Example:

$MDL Klein GDP equation
EQ> y = cn + i + g - t

$STOCHSIMULATE argument
StochStructure = 
list(
  y = list(
      TSRANGE=TRUE,
      TYPE='NORM',
      PARS=c(0,1))
  )
  

then, during the stochastic simulation, the following assignment will be evaluated in the whole simulation TSRANGE:

y <- ( ConstantAdjustment$y + rnorm(...,0,1) ) + cn + i + g - t

EXOGENOUS REFERENCE.

If the generic name var_nameI refers to an exogenous variable, and TYPE!='MATRIX', the disturbance will be additive and applied to the referred exogenous time series.

Example:

$MDL Klein GDP equation
EQ> y = cn + i + g - t

$STOCHSIMULATE argument
StochStructure = 
list(
  g = list(
      TSRANGE=TRUE,
      TYPE='UNIF',
      PARS=c(-1,1))
  )
  

then, during the stochastic simulation, the following assignment will be evaluated in the whole simulation TSRANGE:

y <- ConstantAdjustment$y + cn + i + ( g + runif(...,-1,1) ) - t

If the generic name var_nameI refers to an exogenous variable, and TYPE=='MATRIX', the user-defined disturbance will replace the referred exogenous time series values.

Example:

$MDL equation
EQ> y = cn + i + g - t

$STOCHSIMULATE argument
StochStructure = 
list(
  g = list(
      TSRANGE=myTSRANGE,
      TYPE='MATRIX',
      PARS=userMatrix)
  )
  

then, during the stochastic simulation, the following assignment will be evaluated only in the sub-range myTSRANGE:

y <- ConstantAdjustment$y + cn + i + ( userMatrix ) - t

with userMatrix as a [ myTSRANGE x StochReplica ] matrix.

Generally, each StochStructure list element must be a named list built with the following three element:

- TSRANGE: the time range wherein the disturbance is applied to the referred model variable. 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 disturbance to whole STOCHSIMULATE TSRANGE.

- TYPE: the type of disturbance distribution, .e.g. normal TYPE='NORM', or uniform TYPE='UNIF', or user-defined TYPE='MATRIX';

- PARS: the parameters that shape the disturbance.
In the case of a TYPE='NORM' distribution, these parameters must contain normal distribution mean and standard deviation, e.g. PARS=c(mean,sd);
in the case of a TYPE='UNIF' distribution, these parameters must contain uniform distribution lower and upper bound, e.g. PARS=c(min,max);
in the case of a TYPE='MATRIX', this parameter must contain the user-defined matrix of realized disturbances, e.g. PARS=matrix, with matrix as a [ TSRANGE x StochReplica ] matrix. As said above, if the list element name, e.g. var_nameI, refers to an exogenous variable, the user-built provided matrix will replace the exogenous values in the specified TSRANGE; otherwise, if the element name refers to an endogenous variable, matrix data will be additive to the referred add-factor.

See examples in order to learn how to build a compliant stochastic structure.

StochReplica

an integer value that sets the number of stochastic realizations to be produced

StochSeed

a number used to initialize the pseudo-random number generator. It can be helpful in order to replicate stochastic results

quietly

see SIMULATE

RESCHECKeqList

see SIMULATE

JACOBIAN_SHOCK

see SIMULATE

JacobianDrop

see SIMULATE

forceForwardLooking

see SIMULATE

avoidCompliance

see SIMULATE

...

see SIMULATE

Value

This function will add, into the output BIMETS model object, three new named elements, respectively stochastic_simulation, simulation_MM and INSTRUMENT_MM.

The stochastic_simulation element is a named list(), having endogenous variables as names. Each element will contain two time series: the mean and the standard deviation of the related stochastic simulated endogenous time series.

The arguments passed to the function call during the latest STOCHSIMULATE run will be inserted into the '__STOCH_SIM_PARAMETERS__' element of the stochastic_simulation list; this data can be helpful in order to replicate the stochastic simulation 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 simulation 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 (see example).

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 simulation 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
SIMULATE
RENORM
OPTIMIZE
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 
COMMENT> 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 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 disturbances to exogenous Government Expenditure 'g'
#in whole TSRANGE
myStochStructure <- list(
  cn=list(
        TSRANGE=c(1942,1,1942,1),
        TYPE='NORM',
        PARS=c(0,advancedKleinModel$behaviorals$cn$statistics$StandardErrorRegression)
        ),
  g=list(
        TSRANGE=TRUE,
        TYPE='UNIF',
        PARS=c(-1,1)
        )
  )

#we need to extend exogenous variables up to 1944
advancedKleinModel$modelData <- within(advancedKleinModel$modelData,{
    w2    = TSEXTEND(w2,  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')
  })

#stochastic model forecast
advancedKleinModel <- STOCHSIMULATE(advancedKleinModel
                      ,simType='FORECAST'
                      ,TSRANGE=c(1941,1,1944,1)
                      ,StochStructure=myStochStructure
                      ,StochSeed=123
                      )
                      
#print mean and standard deviation of forecasted GNP
with(advancedKleinModel$stochastic_simulation,TABIT(y$mean, y$sd))

#      Date, Prd., y$mean         , y$sd           
#
#      1941, 1   ,  125.5045      ,  4.250935      
#      1942, 1   ,  173.2946      ,  9.2632        
#      1943, 1   ,  185.9602      ,  11.87774      
#      1944, 1   ,  141.0807      ,  11.6973  
      
#print the unperturbed forecasted GNP along with the
#first 5 perturbed realizations
with(advancedKleinModel$simulation_MM,print(y[,1:6]))

####################################################
#EXAMPLE WITH TYPE='MATRIX'

TSRANGE <- c(1935,1,1940,1)
StochReplica <- 100

#we will perturb simulation by using regression residuals
#get cn and i residuals in TSRANGE
cn_residuals <- TSPROJECT(advancedKleinModel$behaviorals$cn$residuals, 
                          TSRANGE=TSRANGE,
                          ARRAY = TRUE)
i_residuals <- TSPROJECT(advancedKleinModel$behaviorals$i$residuals,
                         TSRANGE=TSRANGE,
                         ARRAY = TRUE)

#define stochastic matrices
cn_matrix <- c()
i_matrix <- c()

#populate matrices
for (idx in 1:StochReplica)
{
  rand <- rnorm(1,0,1)
  cn_matrix <- cbind(cn_matrix,rand*cn_residuals)
  i_matrix <- cbind(i_matrix,rand*i_residuals)
}

#define stochastic structure
myStochStructure <- list(
  cn=list(
    TSRANGE=TRUE,
    TYPE='MATRIX',
    PARS=cn_matrix
  ),
  i=list(
    TSRANGE=TRUE,
    TYPE='MATRIX',
    PARS=i_matrix
  )
)

#stochastic simulation
advancedKleinModel <- STOCHSIMULATE(advancedKleinModel
                                    ,TSRANGE=TSRANGE
                                    ,StochStructure=myStochStructure
)

#print GNP mean and sd
with(advancedKleinModel$stochastic_simulation,TABIT(y$mean, y$sd))

#########################################################
#EXAMPLE OF MODEL THAT REQUIRES THE FULL NEWTON ALGORITHM
#see profit equation

myFullNewtonDefinition<-
"MODEL 

COMMENT> Consumption
BEHAVIORAL> cn
TSRANGE 1922 1 1929 1
EQ> cn =  a1 + a2*p + a3*LAG(p,1) + a4*(w1+w2+w3) 
COEFF> a1 a2 a3 a4

COMMENT> Investment
BEHAVIORAL> i
TSRANGE 1922 1 1929 1
EQ> i = b1 + b2*p + b3*LAG(p,1) + b4*LAG(k,1)
COEFF> b1 b2 b3 b4

COMMENT> Demand for Labor
BEHAVIORAL> w1 
TSRANGE 1922 1 1929 1
EQ> w1 = c1 + c2*(z+y+t-w2) + c3*LAG(z+y+t-w2,1) + c4*time
COEFF> c1 c2 c3 c4

COMMENT> Demand for Labor
BEHAVIORAL> w3 
TSRANGE 1922 1 1929 1
EQ> w3 = c1 + c2*(z+y+t-w2) + c3*LAG(z+y+t-w2,1) + c4*time
COEFF> c1 c2 c3 c4

COMMENT> Gross National Product
IDENTITY> y
EQ> y = cn + i + g - t

IDENTITY> z
EQ> z = cn + i + g - t

COMMENT> Profits with cubic dependence on cn
IDENTITY> p
EQ> p = cn^3/1000 + z + y - (w1+w2+w3)

COMMENT> Capital Stock
IDENTITY> k
EQ> k = LAG(k,1) + i

END
"


#define model data
myModelData<-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)
 )
 

#add data to model
myModelData$z <- myModelData$y
myModelData$w3 <- (myModelData$w1)

myFullNewtonModel <- LOAD_MODEL(modelText=myFullNewtonDefinition)
myFullNewtonModel <- LOAD_MODEL_DATA(myFullNewtonModel,myModelData)

myFullNewtonModel <- ESTIMATE(myFullNewtonModel)
 
#simple Newton will fail, due to 
#large variance in normal disturbances                       
#...while full Newton will converge
myFullNewtonModel <- STOCHSIMULATE(myFullNewtonModel,
                       simAlgo='FULLNEWTON',
                       TSRANGE = c(1921, 1, 1923, 1),
                       simConvergence = 1e-5,
                       simIterLimit = 250,
                       StochReplica = 100,
                       StochSeed=123,
                       StochStructure = list(
                         cn=list(
                           TSRANGE=TRUE,
                           TYPE='NORM',
                           PARS=c(0,20)
                         )
                       )
                       ,verbose=TRUE
                       )