Nonlinear Three-Stage Least Squares Regression

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: NLIN3SLS                                            */
 /*   TITLE: Nonlinear Three-Stage Least Squares Regression      */
 /* PRODUCT: IML                                                 */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: ECON REGR                                           */
 /*   PROCS: PRINT IML                                           */
 /*    DATA:                                                     */
 /*                                                              */
 /* SUPPORT: RHD                         UPDATE:                 */
 /*     REF:                                                     */
 /*    MISC: CONVERTED FROM MATRIX TO IML USING MATIML           */
 /****************************************************************/

*-------THREE STAGE NONLINEAR LEAST SQUARES------------------JPS;

*---GENERATE DATA FOR THE MODEL: -----------*
*     Y1 = A1+EXP(A2*Y2+A3*X1)              *
*     Y2 = B1+B2*X2+B3*Y1*Y1                *
* THIS PROGRAM COULD BE ADAPTED FOR OTHER   *
* SYSTEMS. THE LINES MARKED "***" MUST BE   *
* CHANGED FOR THE MODELS BEING USED.        *
*-------------------------------------------*;

DATA A;
      KEEP X1 X2 Y1 Y2 R1 R2;
L:    N+1; X1=RANUNI(123871); X2=RANUNI(123871);ITER=0;Y1=0;Y2=0;
      IF N<500;
*-----SOLVE BY GAUSS-SEIDEL METHOD-----;
LOOP: Y1LAG=Y1;  Y1=1+EXP(-Y2+X1);
      Y2LAG=Y2;  Y2=1+X2+Y1*Y1;
      ITER+1;  R1=Y1-Y1LAG; R2=Y2-Y2LAG;
      IF ITER<20 & (R1=. | R2=.) THEN GOTO LOOP;
      IF ITER<20 & (ABS(R1)>.0001 | ABS(R2)>.0001) THEN GOTO LOOP;
      IF Y1=. | Y2=. THEN GOTO L;
      Y1=Y1+RANNOR(123831)*.01; Y2=Y2+RANNOR(123831)*.01;
      OUTPUT;
      NOUT+1; IF NOUT<25 THEN GOTO L;
PROC PRINT;

PROC IML;

START MAIN;
  *-----GET THE DATA IN-----;
  USE A ;
  READ ALL INTO DATA ;
  X1= DATA[,1];
  X2= DATA[,2];
  Y1= DATA[,3];
  Y2= DATA[,4];
  FREE DATA;
  N= NROW(X1);
  Z= J(N,1,1)||X1||X2;
  ZI= INV(Z`*Z);
  *-----SET INITIAL PARMETER VALUES-----;
  A={ .9 -1 .9}`;
  B={ .9 1.1  .9}`;
  *-----TWO STAGE NONLINEAR LEAST SQUARES-----*
        *-----FIRST EQUATION-----;
  ITER=0;
  LASTSSE=1E20;
  EPS=1;
  GOTO I1;
L1: ITER=ITER+1;
  J1=-( J(N,1) || Y2#T || X1#T );
    *--EVALUATE THE JACOBIAN--;
  ZJ=Z`*J1;
  D1= SOLVE(ZJ`*ZI*ZJ, ZJ`*ZI*Z`*F1);
 *---SOLVE FOR THE CORRECTION--;
  TRY=0;
  LASTSSE=SSE;
  EPS= MAX( ABS(D1));
  OLDA=A;
H1: A=OLDA-D1;
                         *---APPLY THE CORRECTION--;
I1: A1= A[1,];
  A2= A[2,];
  A3= A[3,];
      *---SEPARATE THE PARAMETERS--;
  T= EXP(A2#Y2+A3#X1);
  F1=Y1-A1-T;
                        *--EVALUATE THE FUNCTION--;
  SSE=(F1`*Z)*ZI*(Z`*F1);
            *--CRITERION TO MINIMIZE--;
  IF ( SSE>LASTSSE) THEN  DO;
            *--HALVING LOGIC--;
     TRY=TRY+1;
     IF ( TRY>10) THEN LINK HALF;
     D1=.5#D1;
     GOTO H1;
     END;
  IF ( EPS>.0001 & ITER<12) THEN GOTO L1;
  MSE=F1`*F1/(N-3);
  COVA= INV(ZJ`*ZI*ZJ)#MSE;
  S=DIAG(1/ SQRT( VECDIAG(COVA)));
  CORRA=S*COVA*S;
  PRINT A[L="MODEL 1 PARAMETER ESTIMATES"];
  PRINT MSE[L="MODEL 1 ESTIMATE OF VARIANCE"];
  PRINT COVA[L="MODEL 1 COVARIANCE OF PARAMETER ESTIMATES"];
  PRINT CORRA[L="MODEL 1 CORRELATION OF PARAMETER ESTIMATES"];
  PRINT F1[L="MODEL 1 RESIDUALS"];
  PRINT ITER, D1;
         *-----SECOND EQUATION;
  ITER=0;
  EPS=1;
  LASTSSE=1E30;
  GOTO I2;
L2: ITER=ITER+1;
  J2=-( J(N,1) || X2 || Y1#Y1 );
     *--EVALUATE THE JACOBIAN--;
  ZJ=Z`*J2;
  D2= SOLVE(ZJ`*ZI*ZJ, ZJ`*ZI*Z`*F2);
 *--SOLVE FOR THE CORRECTION--;
  TRY=0;
  LASTSSE=SSE;
  EPS= MAX( ABS(D2));
  OLDB=B;
H2: B=OLDB-D2;
                         *--APPLY THE CORRECTION--;
I2: B1= B[1,];
  B2= B[2,];
  B3= B[3,];
      *--SEPARATE THE PARAMETERS--;
  F2=Y2-B1-B2#X2-B3#Y1#Y1;
           *--EVALUATE THE FUNCTION;
  SSE=(F2`*Z)*ZI*(Z`*F2);
            *--CRITERION TO MINIMIZE--;
  IF ( SSE>LASTSSE) THEN  DO;
            *--HALVING LOGIC--;
     TRY=TRY+1;
     IF ( TRY>10) THEN LINK HALF;
     D2=.5#D2;
     GOTO H2;
     END;
  IF ( EPS>.0001 & ITER<12) THEN GOTO L2;
  MSE=F2`*F2/(N-3);
  COVB= INV(ZJ`*ZI*ZJ)#MSE;
  S=DIAG(1/ SQRT( VECDIAG(COVB)));
  CORRB=S*COVB*S;
  PRINT B[L="MODEL 2 PARAMETER ESTIMATES"];
  PRINT MSE[L="MODEL 2 ESTIMATE OF VARIANCE"];
  PRINT COVB[L="MODEL 2 COVARIANCE OF PARAMETER ESTIMATES"];
  PRINT CORRB[L="MODEL 2 CORRELATION OF PARAMETER ESTIMATES"];
  PRINT F2[L="MODEL 2 RESIDUALS"];
  PRINT ITER, D2;
  *------THIRD STAGE------;
        *--PROCESS COVARIANCE ACROSS EQUATIONS--;
  NEQ=2;
  *--NUMBER OF EQUATIONS--;
  S=(F1||F2)`*(F1||F2)/N;
  PRINT S[L="COVARIANCE ACROSS MODELS"];
  SC=DIAG(1/ SQRT( VECDIAG(S)));
  SCORR=SC*S*SC;
  PRINT SCORR[L="CORRELATION ACROSS MODELS"];
  SI= INV(S);
  PRINT SI[L="INVERSE OF COVARIANCE ACROSS MODELS"];
        *--USE N2SLS ESTIMATES AS STARTING VALUES--;
  B=A//B;
  ITER=0;
  EPS=1;
  LASTSSE=1E30;
  GOTO I3;
        *--SOLVE THE SYSTEM USING GAUSS-NEWTON--;
L3: ITER=ITER+1;
  J=- BLOCK( J(N,1) || Y2#T || X1#T , J(N,1) || X2 || Y1#Y1);
 *--EVALUATE THE JACOBIAN--;
  ZJ=( I(NEQ)@Z`)*J;
  XPX=ZJ`*(SI@ZI)*ZJ;
  D= SOLVE(XPX,ZJ`*(SI@ZI)*( I(NEQ)@Z`)*F);
  TRY=0;
  LASTSSE=SSE;
  EPS= MAX( ABS(D));
  OLDB=B;
H3: B=OLDB-D;
I3: B1= B[1,];
  B2= B[2,];
  B3= B[3,];
  B4= B[4,];
  B5= B[5,];
  B6= B[6,];
                                *--SEPARATE THE PARAMETERS--;
  T= EXP(B2#Y2+B3#X1);
  F = Y1 - B1-T  // Y2 - B4-B5#X2-B6#Y1#Y1;
     *--EVALUATE THE FUNCTION--;
  ZF=( I(NEQ)@Z`)*F;
  SSE=ZF`*(SI@ZI)*ZF;
              *--CRITERION TO MINIMIZE--;
  IF ( SSE>LASTSSE) THEN  DO;
          *--HALVING LOGIC--;
     TRY=TRY+1;
     IF ( TRY>10) THEN LINK HALF;
     D=.5#D;
     GOTO H3;
     END;
  IF ( EPS>.0001 & ITER<12) THEN GOTO L3;
  PRINT ITER, D;
  PRINT B[L="N3SLS PARAMETER ESTIMATES"];
  COVB= INV(XPX);
  PRINT COVB[L="N3SLS COVARIANCE OF ESTIMATES"];
  SC=DIAG(1/ SQRT( VECDIAG(COVB)));
  CORRB=SC*COVB*SC;
  PRINT CORRB[L="N3SLS CORRELATION OF ESTIMATES"];
  F= SHAPE(F,0,N)`;
  PRINT F[L="N3SLS RESIDUALS"];
  RETURN;
HALF: PRINT
  " NO IMPROVEMENT IN SSE AFTER 10 HALVINGS, CONVERGENCE ASSUMED";
  EPS=0;
  LASTSSE=1E10;
  RETURN;

FINISH MAIN;

RUN MAIN;