| SIMULATE {bimets} | R Documentation |
Simulation of a BIMETS model
Description
The simulation of an econometric model basically consists in solving the system of the equations describing the model for each time period in the specified time interval. Since the equations may not be linear in the variables, and since the graph derived from the "incidence matrix" may be cyclic, the usual methods based on linear algebra are not applicable. The simulation must be solved by using an iterative algorithm (Users can find the "indicence matrix" definition in the following section "The Optimal Reordering").
BIMETS simulation capabilities support:
- Static simulations: a static multiple equation simulation, in which the historical values for the lagged endogenous variables are used in the solutions of subsequent periods;
- Dynamic simulations: a dynamic simulation, in which the simulated values for the lagged endogenous variables are used in the solutions of subsequent periods;
- Residuals check: a single period, single equation simulation; simulated time series in output are just the computation of the RHS (right-hand-side) of their equation, by using the historical values of the involved time series and by accounting for error autocorrelation and PDLs, if any;
- Forecast simulations: similar to dynamic simulation, but during the initialization of the iterative algorithm the starting values of endogenous variables in a period are set equal to the simulated values of the previous period. This allows the simulation of future endogenous observations, i.e. the forecast;
- Stochastic Simulation: see STOCHSIMULATE;
- Partial or total exogenization of endogenous variables: in the provided time interval (i.e. partial exog.) or in whole simulation time range (i.e. total exog.), the values of the selected endogenous variables can be definitely set equal to their historical values, by excluding their equations from the iterative algorithm of simulation;
- Constant adjustment of endogenous variables (add-factors): adds up a new exogenous time series - the "constant adjustment" - in the equation of the selected endogenous variables.
- Gauss-Seidel and Newton-Raphson simulation algorithms: the Gauss-Seidel algorithm is simple, robust, and works well for many backward-looking macro-econometric models. Equations are evaluated as-is in a proper order until the convergence, if any, is verified on the so called "feedback variables" (Users can find the "feedback variable" definition in the next section "The Optimal Reordering"). It is slower than Newton-Raphson algorithms for a very low convergence criterion, and fails to find a convergence for a small set of econometric models, even when a convergence exists. The Newton-Raphson algorithm allows users to solve a broader set of macro-econometric models than the Gauss-Seidel algorithm. Moreover, it is usually faster than the Gauss-Seidel algorithm (on modern computers, users must simulate an extensive econometric model with a low convergence criterion to appreciate the speedup). This type of algorithm requires the construction and the inversion of the Jacobian matrix for the feedback variables; thus, in some scenarios, numerical issues can arise, and users are required to manually exclude some feedback variables from the Jacobian matrix by using the JacobianDrop argument of the SIMULATE procedure.
In details, the generic model suitable for simulation in BIMETS can be written as:
y_1=f_1(\bar{x},\bar{y})
...
y_n=f_n(\bar{x},\bar{y})
being:
n the number of equations in the model;
\bar{y}=[y_1, ... , y_n] the n-dimensional vector of the endogenous variables;
\bar{x}=[x_1, ... , x_m] the m-dimensional vector of the exogenous variables;
f_i(...), i=1..n any kind of functional expression able to be written by using the MDL syntax;
As described later on, in BIMETS a modified Gauss-Seidel iterative algorithm, or a Newton-Raphson algorithm, can solve the system of equations. The convergence properties may vary depending on the model specifications. In some conditions, the algorithm may not converge for a specific model or a specific set of data.
A convergence criterion and a maximum number of iterations to be performed are provided by default. Users can change these criteria by using the simConvergence and simIterLimit arguments of the SIMULATE function.
The general conceptual scheme of the simulation process (for each time period) is the following:
1. initialize the solution for the current simulation period;
2. iteratively solve the system of equations;
3. save the solution, if any;
Step 2 means that for each iteration, the operations are:
2.1 update the values of the current endogenous variables;
2.2 verify that the convergence criterion is satisfied or that the maximum number of allowed iterations has been reached;
The initial solution for the iterative process (step 1) can be given alternatively by:
- the historical value of the endogenous variables for the current simulation period (the default);
- the simulated value of the endogenous variables from the previous simulation period (this alternative is driven by the simType='FORECAST' argument of the SIMULATE function);
In the "dynamic" simulations (i.e. simulations performed by using either the default simType='DYNAMIC' or the simType='FORECAST'), whenever lagged endogenous variables are needed in the computation, the simulated values of the endogenous variables \bar{y} assessed in the previous time periods are used. In this case, the simulation results in a given time period depend on the simulation results in the previous time periods. This kind of simulation is defined as "multiple equation, multiple period".
As an alternative, the actual historical values can be used in the "static" simulations (i.e. simulations performed by using simType='STATIC') rather than simulated values whenever lagged endogenous variables are needed in the computations. In this case, the simulation results in a given time period do not depend on the simulation results in the previous time periods. This kind of simulation is defined as "multiple equation, single period".
The last simulation type available is the residual check (simType='RESCHECK'). With this option, a "single equation, single period" simulation is performed. In this case, no iteration must be carried out. The endogenous variables are assessed for each time period by using historical values for each variable on the right-hand side of the equation, for both lagged and current periods. This kind of simulation helps debug and check of the logical coherence of the equations and the data, and can be used as a simple tool to compute the add-factors.
The debugging of the logical coherence of equations and data is carried out through a Residual Check procedure.
It consists of the following steps:
1. add another exogenous variable - the constant adjustment - to every equation of the model, both behavioral and technical identity: that can be done in BIMETS by using the ConstantAdjustment argument of the SIMULATE function, as in step 3;
2. fill in with the estimated residuals all the constant adjustments for the behavioral equations, and fill in with zeroes the constant adjustments for the technical identities: that can be done in BIMETS by using the SIMULATE procedure with the option simType='RESCHECK', then by analyzing and using the
ConstantAdjustmentRESCHECK attribute of the simulated model, as in the following simulation in step 3.
3. perform a simulation of the model: that can be done in BIMETS by using the SIMULATE procedure with the option
ConstantAdjustment=ConstantAdjustmentRESCHECK;
4. compute the difference between the historical and the simulated values for all the endogenous variables;
5. check whether all the differences assessed in step 4 are zero in whole time range, eventually accounting for the error autocorrelation in behaviorals.
An example on ConstantAdjustmentRESCHECK usage is available at the end of the SIMULATE help page;
If a perfect tracking of the history is obtained, then the equations have been written coherently with the data, otherwise a simulated equation not tracking the historical values is an unambiguous symptom of data inconsistent with the model definition.
Aside from the residual check, the add-factors constitute an important tool to significantly improve the accuracy of forecasts made through an econometric model. Considering the following model:
y_1=f_1(\bar{x},\bar{y}) + z_1
...
y_n=f_n(\bar{x},\bar{y}) + z_n
the add-factors \bar{z}=[z_1, ... ,z_n] can be interpreted as estimates of the disturbance terms' future values or as adjustments of the intercepts in each equation. These add-factors round out the forecasts, by summarizing the effects of all the systematic factors not included in the model. One choice for the computation of the add-factors is given by past estimation residuals and past forecast errors or by an average of these errors.
Please note that, in the case of equation that presents an LHS function, the add-factor will be added before the application of the inverse function, i.e., during the simulation, the following:
g_1(y_1)=f_1(\bar{x},\bar{y}) + z_1
...
g_n(y_n)=f_n(\bar{x},\bar{y}) + z_n
will be solved as:
y_1=g_{1}^{-1}(f_1(\bar{x},\bar{y}) + z_1)
...
y_n=g_{n}^{-1}(f_n(\bar{x},\bar{y}) + z_n)
If a linear dependence between the simulated endogenous and the add-factor is preferred, users can manually insert an auxiliary equation w_i into the model definition, e.g. the following:
g_i(y_i)=f_i(\bar{x},\bar{y})
can be replaced by:
w_i=f_i(\bar{x},\bar{y})
y_i=g_{i}^{-1}(w_i)
During the simulation, the add-factors (if requested by the user) will be applied as in the following:
w_i=f_i(\bar{x},\bar{y}) + v_i
y_i=g_{i}^{-1}(w_i) + z_i
given v_i, z_i as add-factors and the linear dependence from z_i and y_i.
| RATIONAL EXPECTATIONS MODELS |
BIMETS classifies a model as a forward-looking model if any model equation contains the TSLEAD time series function. Forward-looking models assume that economic agents have complete knowledge of an economic system and calculate the future value of economic variables correctly according to that knowledge. Thus, forward-looking models are called also rational expectations models and, in macro-econometric models, model-consistent expectations.
In forward-looking models, simulation complexity arises, and all simulation periods must be solved simultaneously because equations can contain references to past and future values. Thus, given N simulation periods requested by the user, each model equation must be replicated N-1 times and modified before the simulation takes place, accounting for lead transformations. Finally, the extended model must be simulated as a single block of equations.
Internal data structures too, such as the incidence and the Jacobian matrix, and the reordering arrays vpre and vblocks (described later in the "The Optimal Reordering" section), grow with the number of simulation periods requested by the user. Therefore, they can only be calculated when a new simulation is requested rather than when the model MDL definition is parsed, further extending computational time in simulation.
A simulation that spans several decades in a forward-looking models with hundreds of equations is not feasible in BIMETS. For a real scenario in a rational expectations model, see "Computational details and capabilities" section in the "US Federal Reserve quarterly model (FRB/US) in R with bimets" vignette.
To understand BIMETS internals when dealing with forward-looking models, please consider the following simple example of a forward-looking model having a single identity:
IDENTITY> Y EQ> Y = TSLEAD(Y) - TSLAG(Y) + X
Given X as an exogenous variable, if the requested simulation has a TSRANGE that spans two periods, then the model will be internally transformed into something like:
IDENTITY> Y EQ> Y = Y__LEAD__1 - TSLAG(Y) + X IDENTITY> Y__LEAD__1 EQ> Y__LEAD__1 = TSLEAD(Y,2) - Y + TSLEAD(X)
Accordingly, the model will be simulated only on the first period of the TSRANGE. Please note that TSLAG(Y) in the first equation, and TSLEAD(Y,2) in the second equation, are a kind of exogenous variables and must be defined in order for the simulation to be completed. Moreover, Y and Y__LEAD__1 are simultaneously involved in the iterative simulation algorithm, and both depend on each other, as also stated in the incidence matrix for the extended model:
$incidence_matrix
Y Y__LEAD__1
Y 0 1
Y__LEAD__1 1 0
Due to the mechanism described above, only DYNAMIC simulations are allowed in forward-looking models. See examples below, for a Klein-like forward-looking model simulation.
| THE OPTIMAL REORDERING |
In fact, the simulation process takes advantage of an appropriate equations reordering to increase the performances by iteratively solving only one subset of equations, while the others are solved straightforwardly. "...a different ordering of the equations can substantially affect the speed of convergence of the algorithm; indeed some orderings may produce divergence. The less feedback there is, the better the chances for fast convergence..." - Don, Gallo - Solving large sparse systems of equations in econometric models - Journal of Forecasting 1987.
For backward-looking models, the LOAD_MODEL function builds the model's incidence matrix, then defines the proper equations reordering. The incidence matrix is built from the equations of the model; it is a square matrix in which each row and each column represent an endogenous variable. If the (i,j) element is equal to 1 then in the model definition the current value of the endogenous variable referred by the i-row depends directly from the current value of the endogenous variable referred by the j-column. The reader can see an incidence matrix example in the section "BIMETS package" of this manual wherein the content of the kleinModel$incidence_matrix variable is printed out.
In econometric models, the incidence matrix is usually very sparse. Only a few of the total set of endogenous variables are used in each equation. In this situation, ordering the equation in a particular sequence will lead to a sensible reduction of the number of iterations needed to achieve convergence. Reordering the equations is equivalent to rearranging rows and columns of the incidence matrix. In this way, the incidence matrix might be made lower triangular for a subset of the equations.
For this subset, an endogenous variable determined in a specific equation has no incidence in any equation above it, although the same variable might have incidence in equations below it. Such a subset of equations is called recursive. Recursive systems are easy to solve. It is only necessary to solve each equation once if this is done in the proper order. On the other hand, it is unlikely for whole model to be recursive. Indeed the incidence graph is often cyclic, as in the Klein's model that presents the following circular dependecies in the incidence matrix: p <- w1 <- y <- i <- p as shown in the "BIMETS package" section figure.
For some subsets of the equations, some 1's will occur in the upper triangle of the incidence matrix for all possible orderings. Such subsets of equations are called simultaneous. To solve the endogenous variables in the simultaneous subset of equations, an iterative algorithm has to be used. Nevertheless, the equations in a simultaneous subset may be ordered so that the pattern of the 1's in the upper triangle forms a spike. The variables corresponding to the 1's in the upper triangle are called feedback variables.
A qualitative graphical example of an ordered incidence matrix is given in the following figure. The white areas are all 0's, the gray areas contain 0's and 1's. The 1's in the light gray areas refer to variables already evaluated in previous subset of equations, therefore they are known terms within the current subset. The 1's in the dark gray areas refer to variables evaluated within the subset.
In BIMETS, the final pattern of an incidence matrix after the equations reordering generally features N+1 blocks:
- One recursive subset of equation, i.e. the pre-recursive VPRE in image;
- N blocks of equations, VBLOCK in image, each built with a simultaneous VSIM and a post-recursive VPOST subset of equations;
As said, the pre-recursive and the post-recursive subsets are lower triangular. Therefore the corresponding equations are solvable with a cascade substitution with no iteration. Only the simultaneous subsets need an iterative algorithm to be solved. It is important to say that the convergence criterion may also be applied to feedback variables only: when the feedback variables converge, the rest of the simultaneous variables also do.
BIMETS builds and analyzes the model's incidence matrix, and then it i) computes the strongly connected component of the related incidence graph by using the Tarjan algorithm (Ref: Tarjan, Robert - Depth-first search and linear graph algorithms - SIAM Journal on Computing - June 1972), and ii) orders the equations in pre-recursive and, for each block of equations, in simultaneous and post-recursive subsets. The simultaneous subsets are then analyzed in order to find a minimal set of feedback variables. This last problem is known to be NP-complete (Ref: Garey, Johnson - Computers and Intractability: a Guide to the Theory of NP-completeness - San Francisco, Freeman 1979).
The optimal reordering of the model equations is programmatically achieved through the use of an iterative algorithm applied to the incidence matrix that can produce 1+3*N ordered lists of endogenous variables, respectively:
1. One list vpre that is the ordered list containing the names of the endogenous pre-recursive variables to be sequentially computed (once per simulation period) before the simulation iterative algorithm takes place;
2. For each of the N elements in the vblocks list:
2.1 One list vsim (the simultaneous subset) that is the ordered list containing the names of the endogenous variables to be sequentially computed during each iteration of the simulation iterative algorithm;
2.2 One list vfeed that is the list containing the names of the endogenous feedback variables; generally vfeed are the last variables in the ordered vsim list in the sambe block;
2.3. One list vpost that is the ordered list containing the names of the endogenous post-recursive variables to be sequentially computed (once per simulation period) after the simulation iterative algorithm has found a solution in the previous simultaneous subset in the same block;
Once equations are reordered, the previous conceptual scheme is modified as follows:
1. initialize the solution for the current simulation period;
2. compute the pre-recursive equations (i.e. the equation of the endogenous variables in the vpre ordered list);
For each block in vblocks:
3.1 iteratively compute the system of simultaneous equations (i.e. the equation of the endogenous variables in the vsim ordered list): for each iteration: i) update the values of the current endogenous variables, ii) update the feedback variables accordingly to the simulation algorithm in use (see next section for details on simulation algorithms) and iii) verify that the convergence criterion is satisfied on the feedback variables vfeed or that the maximum number of iterations has been reached;
3.2 compute the post-recursive equations (i.e. the equation of the endogenous variables in the vpost ordered list);
Finally:
4. save the solutions;
Clearly, each endogenous variable is updated accordingly to its related equation EQ> in the MDL model definition.
In forward-looking models, the incidence matrix and the equations reordering depend on the simulation periods count, therefore the model attributes incidence_matrix, vblocks, and vpre are calculated only after a simulation has been initialized, and will be available to users in the model$simulation[['__SIM_PARAMETERS__']] lists.
The reader can see an equations reordering example in the section "BIMETS package" of this manual wherein the content of the kleinModel$vpre and kleinModel$vblocks variables are printed out.
| THE SIMULATION ALGORITHMS |
Given x_{j} the j-exogenous variable, j=1..m, and y_{i,k} the i-endogenous variable in a simultaneous subset, at the iteration k, with i=1..n the position of the equation in a reordered model, the modified Gauss-Seidel method takes for the approximation of the endogenous variable y_{i,k} the solution of the following:
y_{i,k}=f_i(x_1, ..., x_m, y_{1,k}, ..., y_{i-1,k}, y_{i,k-1},..., y_{n,k-1})
Newton-Raphson's methods can be seen as an extension of the modified Gauss-Seidel algorithm, and a further step is required: in Newton-Raphson, feedback variables are updated not by using their model equations, but by using the inverse of the Jacobian matrix and the following assignment:
\bar{y}^F_{k} \leftarrow \bar{y}^F_{k-1}+(I-J)^{-1}[\bar{y}^F_{k}-\bar{y}^F_{k-1}]
given the vector of feedback variables values \bar{y}^F_{k}=[y_{n-F+1,k},...,y_{n,k}] at iteration k, the identity matrix I, and the Jacobian matrix J, with I,J \in R^{F,F} and F equal to the number of feedback variables for the given block of equations. Please note that the modified Gauss-Seidel algorithm can be seen as a reduced form of a Netwotn algorithm, given J=0.
In Newton-Raphson methods, the Jacobian matrix J is calculated as follows:
1 - shock the feedback variables one at a time by a small amount;
2 - for each shocked feedback variable, evaluate the shocked solution of the simultaneous subset in the current block;
3 - calculate the derivatives (i.e. the column in the Jacobian matrix related to the shocked feedback variable) using the difference quotients between the shocked and the base solution of the simultaneous subset.
As said, the convergence is always tested at each iteration's end on the feedback variables.
Newton-Raphson methods on a reordered model require the calculation of the Jacobian matrix on the feedback endogenous variables, i.e. at least F+2 iterations per simulation period, with F as the number of feedback variables. For large models (i.e. more than 30 feedback variables) if the overall required convergence is greater than 10^{-6} \% the speedup over the Gauss-Seidel method is small or negative, if the Jacobian matrix is recalculated at each iteration. Moreover, the Gauss-Seidel method does not require a matrix inversion; therefore, it is more robust against algebraical and numerical issues. For small models, both methods are fast on modern computers. On the other hand, Gausse-Seidel fails to find a convergence for a small set of econometric models, even when a convergence exists. In general, given a system of equations Ax=b, with x,b \in R^n, n>0 and A \in R^{n, n}, the Gauss-Seidel algorithm is known to converge if either:
- A is symmetric positive-definite;
- A is strictly or irreducibly diagonally dominant.
To improve simulation speed, BIMETS evaluates the Newton-Raphson algorithm's performance during simulation, and, at each iteration, a new Jacobian matrix is calculated only if the convergence speed is slower than a predefined threshold. In a vectorized simulation (e.g., STOCHSIMULATE, OPTIMIZE, RENORM), if simAlgo="NEWTON" the Jacobian matrix is calculated only on the unperturbed model, then applied to all realizations; if simAlgo="FULLNEWTON" a new Jacobian matrix is calculated for each realization.
The simulation of a non-trivial model, if computed by using the same data but on different hardware, software or numerical libraries, produces numerical differences. Therefore a convergence criterion smaller than 10^{-7} \% frequently leads to a local solution.
See Numerical methods for simulation and optimal control of large-scale macroeconomic models - Gabay, Nepomiastchy, Rachidi, Ravelli - 1980 for further information.
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
SIMULATE( 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,
MULTMATRIX=FALSE,
RENORM=FALSE,
TARGET=NULL,
INSTRUMENT=NULL,
MM_SHOCK=0.00001,
STOCHSIMULATE=FALSE,
StochStructure=NULL,
StochReplica=100,
StochSeed=NULL,
OPTIMIZE=FALSE,
OptimizeBounds=NULL,
OptimizeRestrictions=NULL,
OptimizeFunctions=NULL,
quietly=FALSE,
RESCHECKeqList=NULL,
JACOBIAN_SHOCK=1e-4,
JacobianDrop=NULL,
forceForwardLooking=FALSE,
avoidCompliance=FALSE,
...)
Arguments
model |
The BIMETS model object to be simulated. The simulation requires that all the model behaviorals, if any, have been previously estimated: all the behavioral coefficients (i.e. the regression coefficients and the autoregression coefficients for the errors, if any) must be numerically defined in the model object. (see also |
simAlgo |
The simulation algorithm to be used to solve the system of model equations for each time period in the simulation |
TSRANGE |
The time range of the simulation, as a four dimensional numerical array, |
simType |
The simulation type requested: |
simConvergence |
The percentage convergence value requested for the iterative process, which stops when the percentage difference of all the feedback variables between iterations is less than |
simIterLimit |
The value representing the maximum number of iterations to be performed. The iterative procedure will stop when |
ZeroErrorAC |
If |
BackFill |
Defined as an |
Exogenize |
A named list that specifies the endogenous variables to be exogenized. During the simulation and inside the provided time range, the exogenized endogenous variables will be assigned to their historical values. List names must be the names of the endogenous variables to be exogenized; each element of this list contains the time range of the exogenization for the related endogenous variable, in the form of a 4-dimensional integer array, i.e. start_year, start_period, end_year, end_period. A list element can also be assigned |
ConstantAdjustment |
A named list that specifies the constant adjustments (i.e. add-factors) to be added to the selected equations of the model. Each constant adjustment can be see as a new exogenous variable added to the equation of the specified endogenous variable. The list names are the names of the involved endogenous variables; each element of this is list contains the time series to be added to the equation of the related endogenous variable. Each provided time series must verify the compliance control check defined in |
verbose |
If |
verboseSincePeriod |
An integer that activates the verbose output, during the iterative process, only after the provided number of simulation periods |
verboseVars |
A |
MULTMATRIX |
It is |
RENORM |
It is |
TARGET |
see |
INSTRUMENT |
see |
MM_SHOCK |
see |
STOCHSIMULATE |
It is |
StochStructure |
The |
StochReplica |
An integer value that sets the number of stochastic simulation replications to be performed. See |
StochSeed |
A number used to initialize the pseudo-random number generator. It can be useful in order to replicate stochastic results. See |
OPTIMIZE |
It is |
OptimizeBounds |
see |
OptimizeRestrictions |
see |
OptimizeFunctions |
see |
quietly |
If |
RESCHECKeqList |
If |
JACOBIAN_SHOCK |
The value of the shock added to feedback variables in the derivative calculation of the Jacobian matrix. The default value is |
JacobianDrop |
The array built with feedback variables names to be excluded from the Jacobian matrix calulation |
forceForwardLooking |
If |
avoidCompliance |
If |
... |
Backward compatibility |
Value
This function will add a new named element simulation into the output BIMETS model object.
The new simulation element is a named list; the names of the simulation list are the names of the endogenous variables of the model; each element of the simulation list contains the simulated time series of the related endogenous variable (see example).
The simulation list also contains the '__SIM_PARAMETERS__' element that contains the arguments passed to the function call during the latest SIMULATE run, e.g. TSRANGE, symType, simConvergence, symIterLimit, Exogenize, etc.: this data can be helpful in order to replicate the simulation results.
In case of a simType='RESCHECK' simulation, a new named element ConstantAdjustmentRESCHECK will be added to the output model. This new element is populated with a list of time series that contains, for each endogenous variable, the tracking residuals time series such that, when using this tracking residuals as add-factors in simulation, the related equation will solve to the trajectory given, for that variable, by its historical data (see example).
See Also
MDL
LOAD_MODEL
ESTIMATE
STOCHSIMULATE
MULTMATRIX
RENORM
OPTIMIZE
TIMESERIES
BIMETS indexing
BIMETS configuration
Examples
#define model
myModelDefinition<-
"MODEL
COMMENT> Klein Model 1 of the U.S. Economy
COMMENT> Consumption
BEHAVIORAL> cn
TSRANGE 1921 1 1941 1
EQ> cn = a1 + a2*p + a3*TSLAG(p,1) + a4*(w1+w2)
COEFF> a1 a2 a3 a4
COMMENT> Investment
BEHAVIORAL> i
TSRANGE 1921 1 1941 1
EQ> i = b1 + b2*p + b3*TSLAG(p,1) + b4*TSLAG(k,1)
COEFF> b1 b2 b3 b4
COMMENT> Demand for Labor
BEHAVIORAL> w1
TSRANGE 1921 1 1941 1
EQ> w1 = c1 + c2*(y+t-w2) + c3*TSLAG(y+t-w2,1) + c4*time
COEFF> c1 c2 c3 c4
COMMENT> Gross National Product
IDENTITY> y
EQ> y = cn + i + g - t
COMMENT> Profits
IDENTITY> p
EQ> p = y - (w1+w2)
COMMENT> Capital Stock
IDENTITY> k
EQ> k = TSLAG(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)
)
#load model and model data
myModel<-LOAD_MODEL(modelText=myModelDefinition)
myModel<-LOAD_MODEL_DATA(myModel,myModelData)
#estimate model
myModel<-ESTIMATE(myModel, quietly = TRUE)
#DYNAMIC SIMULATION
#simulate model
myModel<-SIMULATE(myModel
,TSRANGE=c(1923,1,1941,1)
,simConvergence=0.00001
,simIterLimit=100
)
#
#Simulation: 100.00%
#...SIMULATE OK
#get simulated time series "cn" and "y"
TABIT(myModel$simulation$cn)
#
# Date, Prd., myModel$simulation$cn
#
# 1923, 1 , 50.338
# 1924, 1 , 55.6994
# 1925, 1 , 56.7111
# ...
# 1940, 1 , 66.7799
# 1941, 1 , 75.451
#
TABIT(myModel$simulation$y)
#
# Date, Prd., myModel$simulation$y
#
# 1923, 1 , 56.0305
# 1924, 1 , 65.8526
# 1925, 1 , 64.265
# ...
# 1940, 1 , 76.8049
# 1941, 1 , 93.4459
#
#get latest simulation parameters
print(myModel$simulation$'__SIM_PARAMETERS__')
#$TSRANGE
#[1] 1923 1 1941 1
#
#$simType
#[1] "DYNAMIC"
#
#$simConvergence
#[1] 1e-05
#
#$simIterLimit
#[1] 100
#
#$ZeroErrorAC
#[1] FALSE
#
#...etc etc
########################################################
#RESCHECK SIMULATION
#simulate model
myModel<-SIMULATE(myModel
,simType='RESCHECK'
,TSRANGE=c(1923,1,1941,1)
,simConvergence=0.00001
,simIterLimit=100
)
#
#Simulation: 100.00%
#...SIMULATE OK
#get consumption simulated vs historical differences
TABIT(myModel$simulation$cn-myModel$modelData$cn)
#
# Date, Prd., myModel$simulation$cn - myModel$modelData$cn
#
# 1923, 1 , 1.56574
# 1924, 1 , 0.493503
# 1925, 1 , -0.0076079
# ...
# 1939, 1 , -0.989201
# 1940, 1 , -0.785077
# 1941, 1 , 2.17345
#
########################################################
#FORECAST GNP in 1942 and 1943
#we need to extend exogenous variables in 1942 and 1943
myModel$modelData$w2 <- TSEXTEND(myModel$modelData$w2, UPTO=c(1943,1))
myModel$modelData$t <- TSEXTEND(myModel$modelData$t, UPTO=c(1943,1))
myModel$modelData$g <- TSEXTEND(myModel$modelData$g, UPTO=c(1943,1))
myModel$modelData$time <- TSEXTEND(myModel$modelData$time,UPTO=c(1943,1)
,EXTMODE='LINEAR')
#simulate model
myModel<-SIMULATE(myModel
,simType='FORECAST'
,TSRANGE=c(1940,1,1943,1)
,simConvergence=0.00001
,simIterLimit=100
)
#
#Simulation: 100.00%
#...SIMULATE OK
#get forecasted GNP
TABIT(myModel$simulation$y)
#
# Date, Prd., myModel$simulation$y
#
# 1940, 1 , 74.5781
# 1941, 1 , 94.0153
# 1942, 1 , 133.969
# 1943, 1 , 199.913
#
########################################################
#VERBOSE SIMULATION
myModel<-SIMULATE(myModel
,TSRANGE=c(1923,1,1941,1)
,simConvergence=0.00001
,simIterLimit=100
,verbose=TRUE
,verboseSincePeriod=19
,verboseVars=c('cn')
)
########################################################
#DYNAMIC NEWTON SIMULATION
#WITH EXOGENIZATION AND CONSTANT ADJUSTMENTS
#define exogenization list
#'cn' exogenized in 1923-1925
#'i' exogenized in whole TSRANGE
exogenizeList<-list(
cn = c(1923,1,1925,1),
i = TRUE
)
#define add-factors list
constantAdjList<-list(
cn = TIMESERIES(1,-1,START=c(1923,1),FREQ='A'),
y = TIMESERIES(0.1,-0.1,-0.5,START=c(1926,1),FREQ='A')
)
#simulate model
myModel<-SIMULATE(myModel
,simAlgo='NEWTON'
,simType='DYNAMIC'
,TSRANGE=c(1923,1,1941,1)
,simConvergence=0.00001
,simIterLimit=100
,Exogenize=exogenizeList
,ConstantAdjustment=constantAdjList
)
#SIMULATE(): endogenous variable "cn" has been exogenized from (trunc) ...
#SIMULATE(): endogenous variable "i" has been exogenized from (trunc) ...
#SIMULATE(): endogenous variable "cn" has a constant adjustment from (trunc) ...
#SIMULATE(): endogenous variable "y" has a constant adjustment from (trunc) ...
#
#Simulation: 100.00%
#...SIMULATE OK
########################################################
#EXAMPLE OF MODEL THAT FAILS GAUSS CONVERGENCE
#define model
myNewtonModelDefinition<-
"
MODEL
COMMENT> Consumption
BEHAVIORAL> cn
TSRANGE 1922 1 1931 1
EQ> cn = a1 + a2*p + a3*LAG(p,1) + a4*(w1+w2)
COEFF> a1 a2 a3 a4
COMMENT> Investment
BEHAVIORAL> i
TSRANGE 1922 1 1931 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 1931 1
EQ> w1 = 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
COMMENT> Simple copy of y in z
IDENTITY> z
EQ> z = cn + i + g - t
COMMENT> Profits
IDENTITY> p
EQ> p = z + y - (w1+w2)
COMMENT> Capital Stock
IDENTITY> k
EQ> k = LAG(k,1) + i
END
"
#add data to model
myModelData$z <- myModelData$y
myNewtonModel<-LOAD_MODEL(modelText=myNewtonModelDefinition)
myNewtonModel<-LOAD_MODEL_DATA(myNewtonModel,myModelData)
#estimate model
myNewtonModel<-ESTIMATE(myNewtonModel, quietly = TRUE)
#GAUSS simulation fails to converge...
myNewtonModel <- SIMULATE(myNewtonModel,
TSRANGE = c(1921, 1, 1930, 1),
simConvergence = 1e-7)
#...while NEWTON converges
myNewtonModel <- SIMULATE(myNewtonModel,
simAlgo='NEWTON',
TSRANGE = c(1921, 1, 1930, 1),
simConvergence = 1e-7)
########################################################
#EXAMPLE OF MODEL THAT REQUIRES
#A VARIABLE EXCLUSION FROM JACOBIAN MATRIX
#define model
myNewtonWithDropModelDefinition<-
"
MODEL
COMMENT> Consumption
BEHAVIORAL> cn
TSRANGE 1922 1 1931 1
EQ> cn = a1 + a2*p + a3*LAG(p,1) + a4*(w1+w2)
COEFF> a1 a2 a3 a4
COMMENT> Investment
BEHAVIORAL> i
TSRANGE 1922 1 1931 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 1931 1
EQ> w1 = 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
COMMENT> Simple copy of y in z
IDENTITY> z
EQ> z = cn + i + g - t
COMMENT> Profits
IDENTITY> p
EQ> p = z + y - (w1+w2)
IF> y < 0
COMMENT> Capital Stock
IDENTITY> k
EQ> k = LAG(k,1) + i
END
"
#add data to model
myModelData$z <- myModelData$y
myNewtonModel <- LOAD_MODEL(modelText=myNewtonWithDropModelDefinition)
myNewtonModel <- LOAD_MODEL_DATA(myNewtonModel,myModelData)
#estimate model
myNewtonModel <- ESTIMATE(myNewtonModel, quietly = TRUE)
#"p" variable must be removed from Jacobian because of unverified IF>
myNewtonModel <- SIMULATE(myNewtonModel,
simAlgo='NEWTON',
JacobianDrop='p',
TSRANGE = c(1921, 1, 1930, 1),
simConvergence = 1e-7)
########################################################
#COMPARE FORECAST IN 3 ALTERNATIVE
#EXOGENOUS SCENARIOS
#define model
myModelDefinition <-
"MODEL
COMMENT> Klein Model 1 of the U.S. Economy
COMMENT> Consumption
BEHAVIORAL> cn
TSRANGE 1921 1 1941 1
EQ> cn = a1 + a2*p + a3*TSLAG(p,1) + a4*(w1+w2)
COEFF> a1 a2 a3 a4
COMMENT> Investment
BEHAVIORAL> i
TSRANGE 1921 1 1941 1
EQ> i = b1 + b2*p + b3*TSLAG(p,1) + b4*TSLAG(k,1)
COEFF> b1 b2 b3 b4
COMMENT> Demand for Labor
BEHAVIORAL> w1
TSRANGE 1921 1 1941 1
EQ> w1 = c1 + c2*(y+t-w2) + c3*TSLAG(y+t-w2,1) + c4*time
COEFF> c1 c2 c3 c4
COMMENT> Gross National Product
IDENTITY> y
EQ> y = cn + i + g - t
COMMENT> Profits
IDENTITY> p
EQ> p = y - (w1+w2)
COMMENT> Capital Stock
IDENTITY> k
EQ> k = TSLAG(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)
)
#load model and model data
myModel <- LOAD_MODEL(modelText=myModelDefinition)
myModel <- LOAD_MODEL_DATA(myModel,myModelData)
#estimate model
myModel <- ESTIMATE(myModel, quietly = TRUE)
#create 3 new models for the 3 scenarios
modelScenario1 <- myModel
modelScenario2 <- myModel
modelScenario3 <- myModel
#scenario 1, define exogenous paths
modelScenario1$modelData <- within(modelScenario1$modelData,{
w2 = TSEXTEND(w2, UPTO=c(1943,1))
t = TSEXTEND(t, UPTO=c(1943,1))
g = TSEXTEND(g, UPTO=c(1943,1))
time = TSEXTEND(time,UPTO=c(1943,1)
,EXTMODE='LINEAR')
})
#scenario 2, define exogenous paths
modelScenario2$modelData <- within(modelScenario2$modelData,{
w2 = TSEXTEND(w2, UPTO=c(1943,1))
t = TSEXTEND(t, UPTO=c(1943,1))
g = TSEXTEND(g, UPTO=c(1943,1)
,EXTMODE='LINEAR')
time = TSEXTEND(time,UPTO=c(1943,1)
,EXTMODE='LINEAR')
})
#scenario 3, define exogenous paths
#we also change consumption cn add-factor
modelScenario3$modelData <- within(modelScenario3$modelData,{
w2 = TSEXTEND(w2, UPTO=c(1943,1)
,EXTMODE='MEAN4')
t = TSEXTEND(t, UPTO=c(1943,1))
g = TSEXTEND(g, UPTO=c(1943,1)
,EXTMODE='LINEAR')
time = TSEXTEND(time,UPTO=c(1943,1)
,EXTMODE='LINEAR')
})
constantAdjListScenario3 <- constantAdjList
constantAdjListScenario3$cn[[1941,1]] <- c(1,2,3)
#simulate the 3 models
modelScenario1 <- SIMULATE(modelScenario1
,simAlgo='NEWTON'
,simType='FORECAST'
,TSRANGE=c(1940,1,1943,1)
,simConvergence=1e-5
,simIterLimit=20)
modelScenario2 <- SIMULATE(modelScenario2
,simAlgo='NEWTON'
,simType='FORECAST'
,TSRANGE=c(1940,1,1943,1)
,simConvergence=1e-5
,simIterLimit=20)
modelScenario3 <- SIMULATE(modelScenario3
,simAlgo='NEWTON'
,simType='FORECAST'
,TSRANGE=c(1940,1,1943,1)
,simConvergence=1e-5
,simIterLimit=20
,ConstantAdjustment = constantAdjListScenario3
)
#compare results on GNP
TABIT(modelScenario1$simulation$y,
modelScenario2$simulation$y,
modelScenario3$simulation$y)
########################################################
#EXAMPLE OF MODEL'S TRACKING RESIDUALS INITIALIZATION BY USING
#THE RESCHECK SIMULATION'S OUTPUT VALUE "ConstantAdjusmentRESCHECK"
#define the model with LHS funs
myModel<-'MODEL
COMMENT> Modified Klein Model 1 of the U.S. Economy with PDL,
COMMENT> autocorrelation on errors, restrictions and conditional evaluations
COMMENT> LHS functions on EQ
COMMENT> Exp Consumption
BEHAVIORAL> cn
TSRANGE 1925 1 1941 1
EQ> EXP(cn) = a1 + a2*p + a3*TSLAG(p,1) + a4*(w1+w2)
COEFF> a1 a2 a3 a4
ERROR> AUTO(2)
COMMENT> Log Investment
BEHAVIORAL> i
TSRANGE 1925 1 1941 1
EQ> LOG(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
BEHAVIORAL> w1
TSRANGE 1925 1 1941 1
EQ> w1 = c1 + c2*(TSDELTA(y)+t-w2) + c3*TSLAG(TSDELTA(y)+t-w2,1)+c4*time
COEFF> c1 c2 c3 c4
PDL> c3 1 3
COMMENT> Delta Gross National Product
IDENTITY> y
EQ> TSDELTA(y) = EXP(cn) + LOG(i) + g - t
COMMENT> Profits
IDENTITY> p
EQ> p = TSDELTA(y) - (w1+w2)
COMMENT> Capital Stock with switches
IDENTITY> k
EQ> k = TSLAG(k,1) + LOG(i)
IF> LOG(i) > 0
IDENTITY> k
EQ> k = TSLAG(k,1)
IF> LOG(i) <= 0
END'
#define model data
modelData<-list(
cn=TSERIES(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=TSERIES(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=TSERIES(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=TSERIES(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=TSERIES(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=TSERIES(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=TSERIES(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=TSERIES(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=TSERIES(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=TSERIES(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)
)
#example data transformation
modelData<-within(modelData,{
i=exp(i); #we have LOG(i) in the model MDL definition
cn=log(cn); #we have EXP(cn) in the model MDL definition
y=CUMSUM(y) #we have TSDELTA(y) in the model MDL definition
})
#load model and model data
model<-LOAD_MODEL(modelText=myModel)
model<-LOAD_MODEL_DATA(model,modelData)
#estimate model
model<-ESTIMATE(model)
#get initial tracking residuals in range 1925-1935
#we need to set ZeroErrorAC to TRUE due to error autocorrelation
#in Consumption MDL definition
model<-SIMULATE(model,
TSRANGE=c(1925,1,1935,1),
simType='RESCHECK',
ZeroErrorAC=TRUE)
#get init trac
initTrac<-model$ConstantAdjustmentRESCHECK
#dynamic simulation using initTrac as constant adjustments
#we remove first two periods from simulation range
#due to error autocorrelation in Consumption MDL equation
model<-SIMULATE(model,
TSRANGE=c(1927,1,1935,1),
ConstantAdjustment=initTrac)
#check simulated values are equals to historical ones in simulation range
for (idxV in model$vendog)
{
print(max(abs(model$simulation[[idxV]]-
TSPROJECT(model$modelData[[idxV]],TSRANGE=c(1927,1,1935,1)))))
}
########################################################
#EXAMPLE OF FORWARD-LOOKING KLEIN-LIKE MODEL
#HAVING RATIONAL EXPECTATION ON INVESTMENTS
#define model
kleinLeadModelDefinition<-
"MODEL
COMMENT> Klein Model 1 of the U.S. Economy
COMMENT> Consumption
BEHAVIORAL> cn
TSRANGE 1921 1 1941 1
EQ> cn = a1 + a2*p + a3*TSLAG(p,1) + a4*(w1+w2)
COEFF> a1 a2 a3 a4
COMMENT> Investment with TSLEAD
IDENTITY> i
EQ> i = (MOVAVG(i,2)+TSLEAD(i))/2
COMMENT> Demand for Labor
BEHAVIORAL> w1
TSRANGE 1921 1 1941 1
EQ> w1 = c1 + c2*(y+t-w2) + c3*TSLAG(y+t-w2,1) + c4*time
COEFF> c1 c2 c3 c4
COMMENT> Gross National Product
IDENTITY> y
EQ> y = cn + i + g - t
COMMENT> Profits
IDENTITY> p
EQ> p = y - (w1+w2)
COMMENT> Capital Stock
IDENTITY> k
EQ> k = TSLAG(k,1) + i
END"
#define model data
kleinLeadModelData<-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 model and model data
kleinLeadModel<-LOAD_MODEL(modelText=kleinLeadModelDefinition)
kleinLeadModel<-LOAD_MODEL_DATA(kleinLeadModel,kleinLeadModelData)
#estimate model
kleinLeadModel<-ESTIMATE(kleinLeadModel, quietly = TRUE)
#set expected value of 2 for Investment in 1931
#(note that simulation TSRANGE spans up to 1930)
kleinLeadModel$modelData$i[[1931,1]]<-2
#simulate model
kleinLeadModel<-SIMULATE(kleinLeadModel
,TSRANGE=c(1924,1,1930,1))
#print simulated investments
TABIT(kleinLeadModel$simulation$i)