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;