Single-Equation Solution

For normalized form equation systems, the solution either can simultaneously satisfy all the equations or can be computed for each equation separately, by using the actual values of the solution variables in the current period to compute each predicted value. By default, PROC MODEL computes a simultaneous solution. The SINGLE option in the SOLVE statement selects single-equation solutions.

Single-equation simulations are often made to produce residuals (which estimate the random terms of the stochastic equations) rather than the predicted values themselves. If the input data and range are the same as that used for parameter estimation, a static single-equation simulation reproduces the residuals of the estimation.

Newton’s Method

The NEWTON option in the SOLVE statement requests Newton’s method to simultaneously solve the equations for each observation. Newton’s method is the default solution method. Newton’s method is an iterative scheme that uses the derivatives of the equations with respect to the solution variables, , to compute a change vector as

PROC MODEL builds and solves by using efficient sparse matrix techniques. The solution variables y at the ith iteration are then updated as

where d is a damping factor between 0 and 1 chosen iteratively so that

The number of subiterations allowed for finding a suitable d is controlled by the MAXSUBITER= option. The number of iterations of Newton’s method allowed for each observation is controlled by MAXITER= option. See Ortega and Rheinbolt (1970) for more details.

Jacobi Method

The JACOBI option in the SOLVE statement selects a matrix-free alternative to Newton’s method. This method is the traditional nonlinear Jacobi method found in the literature. The Jacobi method as implemented in PROC MODEL substitutes predicted values for the endogenous variables and iterates until a fixed point is reached. Then necessary derivatives are computed only for the diagonal elements of the jacobian, J.

If the normalized form equation is

the Jacobi iteration has the form

Seidel Method

The Seidel method is an order-dependent alternative to the Jacobi method. The Seidel method is selected by the SEIDEL option in the SOLVE statement. The Seidel method is like the Jacobi method except that in the Seidel method the model is further edited to substitute the predicted values into the solution variables immediately after they are computed. Seidel thus differs from the other methods in that the values of the solution variables are not fixed within an iteration. With the other methods, the order of the equations in the model program makes no difference, but the Seidel method might work much differently when the equations are specified in a different sequence. Note that this fixed point method is the traditional nonlinear Seidel method found in the literature.

The iteration has the form

where is the jth equation variable at the ith iteration and

If the model is recursive, and if the equations are in recursive order, the Seidel method converges at once. If the model is block-recursive, the Seidel method might converge faster if the equations are grouped by block and the blocks are placed in block-recursive order. The BLOCK option can be used to determine the block-recursive form.

Jacobi and Seidel Methods with General Form Equations

Jacobi and Seidel solution methods support general form equations.

There are two cases where derivatives are (automatically) computed. The first case is for equations with the solution variable on the right-hand side and on the left-hand side of the equation

In this case the derivative of ERROR. with respect to is computed, and the new approximation is computed as

The second case is a system of equations that contains one or more EQ. equations. In this case, a heuristic algorithm is used to make the assignment of a unique solution variable to each general form equation. Use the DETAILS option in the SOLVE statement to print a listing of the assigned variables.

Once the assignment is made, the new approximation is computed as

If is the number of general form equations, then derivatives are required.

The convergence properties of the Jacobi and Seidel solution methods remain significantly poorer than the default Newton’s method.

Comparison of Methods

Newton’s method is the default and should work better than the others for most small- to medium-sized models. The Seidel method is always faster than the Jacobi for recursive models with equations in recursive order. For very large models and some highly nonlinear smaller models, the Jacobi or Seidel methods can sometimes be faster. Newton’s method uses more memory than the Jacobi or Seidel methods.

Both the Newton’s method and the Jacobi method are order-invariant in the sense that the order in which equations are specified in the model program has no effect on the operation of the iterative solution process. In order-invariant methods, the values of the solution variables are fixed for the entire execution of the model program. Assignments to model variables are automatically changed to assignments to corresponding equation variables. Only after the model program has completed execution are the results used to compute the new solution values for the next iteration.

Troubleshooting Problems

In solving a simultaneous nonlinear dynamic model you might encounter some of the following problems.

Convergence Problems

The following items can cause convergence problems:

• There are illegal function values ( for example ).

• There are local minima in the model equation.

• No solution exists.

• Multiple solutions exist.

• Initial values are too far from the solution.

• The CONVERGE= value is too small.

When PROC MODEL fails to find a solution to the system, the current iteration information and the program data vector are printed. The simulation halts if actual values are not available for the simulation to proceed. Consider the following program, which produces the output shown in Figure 19.82:

data test1;
   do t=1 to 50;
      x1 = sqrt(t) ;
      y = .;
      output;
   end;

proc model data=test1;
   exogenous x1 ;
   control a1 -1 b1 -29 c1 -4 ;
   y = a1 * sqrt(y) + b1 * x1 * x1 + c1 * lag(x1);
   solve y / out=sim forecast dynamic ;
run;

                              The MODEL Procedure                               
                        Dynamic Single-Equation Forecast                        
                                                                                
                     --- Listing of Program Data Vector ---                     
   _N_:               12     ACTUAL.x1:    1.41421     ACTUAL.y:           .    
   ERROR.y:            .     PRED.y:             .     a1:                -1    
   b1:               -29     c1:                -4     x1:           1.41421    
   y:           -0.00109                                                        
   @PRED.y/@y:           .   @ERROR.y/@y:          .                            
                                                                                

At the first observation, a solution to the following equation is attempted:

There is no solution to this problem. The iterative solution process got as close as it could to making Y negative while still being able to evaluate the model. This problem can be avoided in this case by altering the equation.

In other models, the problem of missing values can be avoided by either altering the data set to provide better starting values for the solution variables or by altering the equations.

You should be aware that, in general, a nonlinear system can have any number of solutions and the solution found might not be the one that you want. When multiple solutions exist, the solution that is found is usually determined by the starting values for the iterations. If the value from the input data set for a solution variable is missing, the starting value for it is taken from the solution of the last period (if nonmissing) or else the solution estimate is started at 0.

proc model data=test1;
   exogenous x1 ;
   control a1 -1 b1 -29 c1 -4 ;
   y = a1 * sqrt(abs(y)) + b1 * x1 * x1 + c1 * lag(x1);
   solve y / out=sim forecast dynamic itprint;
run;