* box cox 10-23-98  - Make sure your data is y x1 x2 x3 and so on.
* If you run this more than once, read in the orginial data
* Make sure that all y >0.

* SAVE YOU DATA BEFORE RUNNING THIS.

SET MXLOOP =9999.


MATRIX.
GET W/VARIABLES = all/FILE=*/missing=omit.

COMPUTE NC = NCOL(W).

COMPUTE NR =NROW(W).

COMPUTE X=MAKE(NR,1,0).
COMPUTE BIGX =MAKE(NR,NC,1).
LOOP II = 1 TO NR.
COMPUTE X(II,1)=W(II,1).
LOOP JJ  = 1 TO NC-1.
COMPUTE BIGX(II,JJ+1)=W(II,JJ+1).
END LOOP.
END LOOP.

* PRINT BIGX.

COMPUTE NL =60.
COMPUTE LAMDA1=MAKE(NL,1,0).
PRINT NC.
PRINT NR.
COMPUTE J=MAKE(NR,1,1).
COMPUTE A = LN(X).
COMPUTE MEANS =T(J)*A/NR.
 *PRINT MEANS.
COMPUTE B=MAKE(NR,NL,0).
* PRINT MEANS.
 COMPUTE MEANS =EXP(MEANS).
*PRINT MEANS.
COMPUTE LAM = -3.1.
LOOP I = 1 TO NL .
COMPUTE LAM = LAM +.1.
COMPUTE LAMDA1(I,1) =LAM.
*PRINT LAM.
DO IF  ABS (LAM) >.001.
COMPUTE A1 = X&**LAM -1.
*PRINT A1.
COMPUTE A2 =(T(LAM*(MEANS)&**(LAM-1))).
* PRINT A2.
COMPUTE Y= A1*INV(A2).
* PRINT Y.
ELSE.
 COMPUTE Y = kroneker(means,LN(X)).
END IF.
LOOP JJ = 1 TO NR.
COMPUTE B(JJ,I)=Y(JJ,1).
END LOOP.
END LOOP.
 * PRINT B . 


COMPUTE BETA = INV(T(BIGX)*BIGX)*T(BIGX)*B.
COMPUTE SS = T(B-BIGX*BETA)*(B-BIGX*BETA)/NR.
COMPUTE RMS =DIAG(SS).
* PRINT RMS.
* PRINT LAMDA1.

COMPUTE ANS=MAKE(NL,2,0).

LOOP I = 1 TO NL.
COMPUTE ANS(I,1) = LAMDA1(I,1).
COMPUTE ANS(I,2)= sqrt(RMS(I,1)).
END LOOP.
* PRINT ANS.
save ANS /outfile=*.
END MATRIX.
rename  variables col1= lambda col2= rmse.
GRAPH
  /SCATTERPLOT(BIVAR)=lambda WITH rmse  /MISSING=LISTWISE
  /TITLE= 'RMSE VERSUS LAMBDA'.