``` /****************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                 */
/*                                                              */
/*    NAME: RSM                                                 */
/*   TITLE: Quadratic Response Surface Regression               */
/* PRODUCT: IML                                                 */
/*  SYSTEM: ALL                                                 */
/*    KEYS: MATRIX  REGR    SUGI6                               */
/*   PROCS: IML                                                 */
/*    DATA:                                                     */
/*                                                              */
/* SUPPORT: RHD                         UPDATE:                 */
/*     REF:                                                     */
/*    MISC:                                                     */
/*                                                              */
/****************************************************************/

proc iml;

| This matrix routine reads in the factor variables and the    |
| response, forms the quadratic regression model and estimates |
| the parameters, then solves for the optimal response, prints |
| the optimal factors and response, and then displays the      |
| the eigenvalues and eigenvectors of the matrix of quadratic  |
| parameter estimates to determine if the solution is a maximum|
| or minimum, or saddlepoint, and which direction has the      |
| steepest and gentlest slopes.                                |
*--------------------------------------------------------------*
| Given that D contains the factor variables,                  |
|        and Y contains the response.                          |
*--------------------------------------------------------------*/
START RSM;
N=NROW(D);  K=NCOL(D);               /* dimensions */
X=J(N,1,1)||D;                       /* set up design matrix */
DO I=1 TO K; DO J=1 TO I; X=X||D[,I]#D[,J]; END; END;

BETA=SOLVE(X`*X,X`*Y);          /* solve parameter estimates */
PRINT "Parameter Estimates" , BETA;

C=BETA[1];                      /* intercept estimate */
B=BETA[2:(K+1)];               /* linear estiamtes */
A=J(K,K,0); L=K+1;             /* form quadratics into matrix */
DO I=1 TO K; DO J=1 TO I; L=L+1; A[I,J]=BETA[L,]; END; END;
A=(A+A`)*.5;                   /* symmetrize */

XX=-.5*SOLVE(A,B);             /* solve for critical value */
PRINT , "Critical Factor Values" , XX;
YOPT=C + B`*XX + XX`*A*XX; /*compute response at critical value*/
PRINT , "Response at Critical Value" YOPT;

CALL EIGEN(EVAL,EVEC,A);
PRINT , "Eigenvalues and Eigenvectors", EVAL, EVEC;

IF MIN(EVAL)>0 THEN PRINT , "Solution Was a Minimum";
IF MAX(EVAL)<0 THEN PRINT , "Solution Was a Maximum";
FINISH;

/* Sample Problem with Two Factors */

D={-1 -1,
-1  0,
-1  1,
0 -1,
0  0,
0  1,
1 -1,
1  0,
1  1};

Y={ 71.7, 75.2, 76.3, 79.2, 81.5, 80.2, 80.1, 79.1, 75.8};

RUN RSM;

```