* box cox 2/28/03  - 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.
* Now has  a " 95% CI for lambda".
* any lambda is acceptable provided the associated rmse <= ci95.

* 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 = gINV(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,3,0).

LOOP I = 1 TO NL.
COMPUTE ANS(I,1) = LAMDA1(I,1).
COMPUTE ANS(I,2)= sqrt(RMS(I,1)).
END LOOP.
compute min1 = cmin(ans).
print min1.
loop i = 1 to nl.
compute ans( i,3)  = min1(2)*exp(3.84/(2*nr)) .
end loop.
* PRINT ANS.
save ANS /outfile=*.

print min1(2).
print ans(1,3).
END MATRIX.
rename  variables col1= lambda col2= rmse  col3 = ci95.
GRAPH
  /SCATTERPLOT(OVERLAY)=lambda lambda  WITH ci95 rmse (PAIR)
  /MISSING=LISTWISE .