* 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

* SAVE YOU DATA BEFORE RUNNING THIS.
new file.
input Program.
compute #n = 2500.
loop i = 1 to #n.
leave i.
end case.
end loop.
end file.
end input program.
execute.

COMPUTE i = rv.beta (2,8) .
EXECUTE .
SAVE OUTFILE='C:\test1.sav'
 /COMPRESSED.


EXAMINE
  VARIABLES=i
  /PLOT BOXPLOT HISTOGRAM NPPLOT
  /COMPARE GROUP
  /STATISTICS DESCRIPTIVES
  /CINTERVAL 95
  /MISSING LISTWISE
  /NOTOTAL.

SET MXLOOP =99999.
set workspace 90000.

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 =200.
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 = -1.1.
LOOP I = 1 TO NL .
COMPUTE LAM = LAM +.01.
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=*.

compute min1 = cmin(ans).
print min1.

COMPUTE D2 = MIN1(1,2).
PRINT D2.

compute lamb1 = ans(1,1).

loop i = 1 to nl.
COMPUTE D1 = ANS(I,2).

 do if   abs (D1-D2) < .0000001 .
compute   lamb1 = ans(i,1).
end if.

end loop.

print lamb1

END MATRIX.
rename  variables col1= lambda col2= rmse.
GRAPH
  /SCATTERPLOT(BIVAR)=lambda WITH rmse  /MISSING=LISTWISE
  /TITLE= 'RMSE VERSUS LAMBDA'.
SAVE OUTFILE='C:\lambda.sav'
 /COMPRESSED.



GET
  FILE='C:\test1.sav'.
COMPUTE LAMBDA = .27.
compute x = i**LAMBDA.
EXAMINE
  VARIABLES=x  /PLOT BOXPLOT HISTOGRAM NPPLOT
  /COMPARE GROUP
  /STATISTICS DESCRIPTIVES
  /CINTERVAL 95
  /MISSING LISTWISE
  /NOTOTAL.


COMPUTE LAMBDA = .27.
COMPUTE mean = 1.68**(1/LAMBDA) .
COMPUTE TM = 10.
COMPUTE bIAS =(100./10.)*LAMBDA + ABS(.25-LAMBDA) .
COMPUTE DIFF = (tM - (MEAN+BIAS)).
EXECUTE .


DESCRIPTIVES
  VARIABLES=i x
  /STATISTICS=MEAN VARIANCE .