SECTION 14

SIMULTANEOUS EQUATIONS ESTIMATION BY FIML (SIMEQ)

SIMEQ permits the estimation of the parameters of both linear and non- linear simultaneous equations systems by full-information maximum likelihood. For linear systems, the equations corresponding to the ith observation are expressed as 

                          Ay(i) + Bx(i) = u(i)           (1)

where y(i) = vector of the ith observation on the jointly dependent variables, x(i) = vector of the ith observation on the predetermined variables, A = a nonsingular matrix of coefficients of order ig x ig, where ig is the number of equations in the system B = a ig x k matrix of coefficients, where k is the number of predetermined variables in the system, u(i) = vector of the ith realization of jointly normal errors. Nonlinear systems are written as 

                         F(y(i),x(i),C) = u(i)           (2)

where F is a vector function and where C represents all unknown parameters. It is assumed that (a) the number of observations is greater than the number of parameters to be estimated, (b) that every equation is identified or overidentified, (c) that all identities have been eliminated.

SECTION 14.1 THE MAIN PROGRAM

The MAIN program must contain the following: 

      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION P(np)
      CHARACTER*8 ALABEL(np)
      COMMON/BSMLQ0/X(n,k),Y(n,ig)
      COMMON/BSMLQ2/IG,NBEG,NEND,IAUTO,ISDIAG,IJAC
      COMMON/BSMLQ3/IQQ2(10,10)
      COMMON/BSMLQ4/MP
C for details about these COMMON blocks see below in this
C subsection and also in Section 14.2
      EXTERNAL SIMEQ1,PERM
C in the next statement "method" refers to GRADX or DFP or etc.
      EXTERNAL method
C statements to read in X and Y
C statements to read in starting values for the elements of P
C statements to read in IQQ2 (if needed)
C statements to set up OPT (set ACC,ITERL, MAX=1,and NPE) (See Sect. 1.1 and 1.11)
      NP=np
C note that NPE .LE. NP
      MP=NP
      IG = ig
      NBEG =
      NEND =
      IAUTO=
      ISDIAG=
      IJAC =
      IP =
      JP =
C Possibly read in IQQ2 array---see below
      CALL LABEL(ALABEL,NP)
      CALL SMLCLN(IER,NP,P,ALABEL,IP,JP)
      CALL OPT(P,NPE,F,METHOD,ITERL,MAX,IER,ACC,SIMEQ1,ALABEL)

where

X = observations on the X-variables (see Eq. (1))
Y = observations on the Y variables (see Eq. (1))
P = vector of parameters to be estimated
NBEG = index of first observation to be used (normally 1)
NEND = index of last observation to be used (normally n)
np is the maximum number of parameters to be estimated,
n is the maximum number of observations on X and Y,
k is the number of actual columns in X,
ig is the number of columns in Y (equal to the number of equations),
number of equations in system; NOTE: IG must not exceed 10 in this version
ISDIAG = 0 if the covariance matrix is not to be diagonal (full and unrestricted)
ISDIAG = 1 if it is to be diagonal
ISDIAG = 2 if selective covariance matrix elements are to be restricted (see IQQ2)
IJAC =1 if the Jacobian array is computed by the user in the SUBROUTINE SIMEQ3 and if the Jacobian is invariant with the observations
IJAC = 2 f the Jacobian array is computed by the users in SIMEQ3 AND the elements of the array are not invariant with observations
IAUTO = 0 if no autocorrelations of residuals are to be estimated
IAUTO = 1 if they are to be estimated
IQQ2 = an ig x ig array; if (ij)th element = 0, the corresponding covariance is restricted to be zero; if = 1, the corresponding covariance elements is unrestricted (see ISDIAG = 2)
JP = 1 causes SMLCLN to call PRMCHK (see Section 1.9)
JP = 0 causes PRMCHK not to be called
IP = 1 (normal setting) is used by PRMCHK to do a permutation automatically

The SUBROUTINE SMLCLN checks for the validity of the input parameters and calls SUBROUTINE PRMCHK if optimization over a subset of the parameters is desired, which is indicated by setting JP = 1 before calling SMLCLN. If JP = 0, PRMCHK will not be called. In that event, IP is of no relevance, since it pertains only to how PRMCHK operates. The normal setting for IP is IP = 1. (See Section 1.9) After returning from SMLCLN the user should check the value of IER and stop computation if an error has been discovered. (See Error Codes at the end of this Section.) If optimization over a subset of parameters is desired, the call to OPT should be 

      CALL OPT(P,NPE,F,METHOD,ITERL,MAX,IER,ACC,PERM,ALABEL)

and the user must also provide a SUBROUTINE PERM as follows:

      SUBROUTINE PERM(A,NPE,FU0,*)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(np)
      EXTERNAL SIMEQ1
      CALL PERM1(A,np,NPE,FU0,&10,SIMEQ1)
      RETURN
10    RETURN 1
      END

SECTION 14.2 DETAILED CONSIDERATIONS AND EXAMPLE

The parameters are stored in an array P. The order of storage is arbitrary with the following exceptions: (a) No storage is provided for the elements of the covariance matrix because the likelihood is maximized in condensed form; (b) Autocorrelation coefficients of the error terms must follow all other coefficients and must be in the order in which the equations are written.

EXAMPLE. Denote the coefficients by c1, c2, etc., and let the equation system be the following:

            y(1) + c2y(2)          = c7x(1) + c8
                     y(2) + c4y(3) = c9x(2) + c10
            y(1)          + c6y(3) = c11x(2) + c12x(3) + c13

It is convenient to consider the coefficient of y(1) in the first equation as the coefficient c1 (= 1.0), the coefficient of y(2) in the second equation as c3 and the coefficient of y(1) in the third equation as c5. Then the P( ) array could be set to represent (c1,c2,c3,c4,c5,c6,c7,c8,c9,c10, c11,c12,c13). If three autocorrelation coefficient were to be estimated, the P( ) array would be (c1,c2,c3,c4,c5,c6,c7,c8,c9,c10, c11,c12,c13,c14,c15,c16), with the last three representing rho1,rho2,rho3.Of course, if any of the coefficients is a constant, the appropriate elements of the P( ) array must be 'PERMED' out, i.e., no optimization must be done with respect to coefficients that are constants. In this case the IPRM array would contain, in sequence, the following values: 2 4 6 7 8 9 10 11 12 13 1 3 5 14 15 16, if no autocorrelations are estimated. In this case NP=16, MP=16, NPE=10. It is particularly important to arrange storage correctly if optimization is done with respect to a subset of the full set of parameters. EXAMPLE. Assume that there is a "full" set of parameters consisting of 16 parameters. Assume that it is a 3-equation system, so that there are 3 autocorrelation coefficients. Further assume that we do not wish to optimize with respect to the 4th, 6th and 11th parameters, but do wish to estimate the autocorrelations. Then, before calling OPT, we should set NP=16 MP=16 NPE=13 and the IPRM array should contain the following numbers in its 16 elements: 1 2 3 5 7 8 9 10 12 13 14 15 16 4 6 11. The autocorrelation coefficients are located in elements 14,15,16. The user must provide the appropriate array for the Jacobian computation by filling in the upper IG x IG portion of the YJAC array in SUBROUTINE SIMEQ3. This can be done once-and-for-all if the elements of the Jacobian do not change with the observations (set IJAC=1) and in this case the parameter I in the calling sequence is irrelevant. For the example at hand the Jacobian array is 

                                  1   c2  0
                                  0   1   c4
                                  1   0   c6

For the arrangement of coefficients in the P( ) array indicated above, the SUBROUTINE SIMEQ3 is as follows: 

      SUBROUTINE SIMEQ3(P,NP,YJAC,IG,I)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION P(NP),YJAC(10,10)
      DO I=1,3
      DO J=1,3
      YJAC(I,J)=0.
      END DO
      END DO
      YJAC(1,1)=1.0
      YJAC(1,2)=P(2)
      YJAC(2,2)=1.0
      YJAC(2,3)=P(4)
      YJAC(3,3)=P(6)
      YJAC(3,1)=1.0
      RETURN
      END

where YJAC is the name of the Jacobian array. In the present example, the input parameter I is irrelevant because the Jacobian does not depend on the observations. If it did, SIMEQ3 would have to be programmed so that it computes the relevant Jacobian for the ith observation. SIMEQ also assumes that the user provides 

      SUBROUTINE SIMEQ2(NP,P,UVEC,I,*)

which computes the residuals of the equations for observation I. For the example above, this is as follows: 

      SUBROUTINE SIMEQ2(NP,P,UVEC,I,*)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION UVEC(10),P(NP)
      COMMON/BSMLQ0/X(60,3),Y(60,3)
      UVEC(1)=Y(I,1)+P(2)*Y(I,2)-P(7)*X(I,1)-P(8)
      UVEC(2)=Y(I,2)+P(4)*Y(I,3)-P(9)*X(I,2)-P(10)
      UVEC(3)=Y(I,1)+P(6)*Y(I,3)-P(12)*X(I,3)-P(11)*X(I,2)-P(13)
      RETURN
      END

If, in executing SIMEQ2, it is possible for errors to occur (e.g., if in a nonlinear model we need to evaluate an expression such as DLOG(P(17)+P(18)*X(I,4)) where the argument of the DLOG function might be negative, the user should test for the error condition and if it does occur, a RETURN 1 should be executed.

SECTION 14.3 ERROR RETURNS

The following error returns are possible: 

          Type of error
IER = -51 IQQ2
IER = -53 IG
IER = -55 IG and IQQ2
IER = -57 (NBEG or NEND)
IER = -59 (NBEG or NEND) and IQQ2
IER = -61 IG and (NBEG or NEND)
IER = -63 IG and (NBEG or NEND) and IQQ2
IER = -65 PRMCHK error in SMLCLN: an invalid permutation was given to PRMCHK.
IER = -66 Jacobian is singular
IER = -67 Covariance matrix is singular
IER = -70 Error return from SIMEQ2
IER = -71 Determinant of covariance matrix is .LE.0
IER = -72 Autocorrelations are not inside (-1.0,1.0)

Errors -66 through -72 are produced in SUBROUTINE SIMEQ1; the function subroutine that it being optimized. Hence these IER values are inaccessible during optimization and may be changed in subsequent iterations. If it is important for the user to know the last error produced prior to exiting from OPT, the user should introduce into this subroutine a labelled COMMON block to pass the value of IER back to the calling program. Note that the existing COMMON/BOPT2/ should not be used for this purpose because it is used by OPT, and OPT may reset the value of IER.

SECTION 14.4 TWO-STAGE LEAST SQUARES (TSLS)

TSLS permits the estimation of a single strcutural equation written as

                          y = Y1γ + X1β + u

where

y is an n x 1 vector of observations on the dependent variable,
Y1 is an n x l1 matrix of observations on other included endogenous variables,
γ is an l1 vector of parameters to be estimated,
X1 is an n x k1 matrix of included exogenous variables,
β is a k1 vector of parameters to be estimated,
u is an n-vector of unobserved error terms.

The call to TSLS is

     CALL TSLS(Y,Y1,X,G,B,ASTD,ZZ,N,L1,K,K1,MEXOG,IER)

where

Y, Y1, X, N, L1, K, K1, MEXOG are inputs, G, B, ASTD, ZZ, IER are outputs, and where Y (y), Y1 (y1, G (γ), B (β), N (n), L1 (l1), K1 (k1) have the meanings given above, and where
ASTD is an L1+K1 vector of asymtptotic standard errors
ZZ is an L1+K1 by L1+K1 array that will contain the estimate of the asymptotic covariance matrix of G and B,
K is the total number of exogeneous variables in the model,
MEXOG is an integer vector of length K containing the identifiers of the included exogenous variables. Thus, if the equation contains on the right hand side the 3rd, 5th and 6th exogenous variables from the array X, we would have MEXOG(1)=3, MEXOG(2)=5, MEXOG(3)=6.

The MAIN program that calls TSLS should contain

      COMMON/BPRINT/IPT,NFILE,NDIG,NPUNCH,JPT,MFILE
      CALL DFLT

Error Returns:

IER 	= 0 normal return
	= -12 error in MEXOG (a MEXOG value of <= 0 or > K)
	= -13 error in MEXOG (two MEXOG values are identical)
	= -16 error in defining allocatable arrays
	= -17 X'X is singular
	= -18 ZZ is singular
	= -19 ZZ is not positive definite.
	
      


Return to the Beginning