** PROGRAM IS BLUP.SAS.  COMPUTES BLUPS FOR SIZE N PROBLEM
PLUS PREDICTION ERROR VARIANCE. **; 
options ls=85;

PROC IML; RESET FUZZ=.000001 fw=5 spaces=5; a=.9; B=.5; 
G=SHAPE(0,1,33); 
G[1,1]=1+(A-B)*(A-B)/(1-A*A); G[1,2]=A*G[1,1]- B ; 
DO LAG1=3 TO 33; G[1,LAG1]=A*G[1,LAG1-1];  END; 

N=32; GAMMA=TOEPLITZ(G[1,1:N]); GIN=INV(GAMMA); 
N1=N+1; VEC=G[1,2:N1]`; BETA=GIN*VEC; 
TAUSQ=G[1,1]-BETA`*GAMMA*BETA; 
COEF=BETA[1:8,1]`//BETA[9:16,1]`//BETA[17:24,1]`//BETA[25:32,1]`; 
PRINT A B COEF TAUSQ; 

DO N=2 to 4;
GAMMA=TOEPLITZ(G[1,1:N]); GIN=INV(GAMMA); 
N1=N+1; VEC=G[1,2:N1]`; BETA=GIN*VEC; 
TAUSQ=G[1,1]-BETA`*GAMMA*BETA; 
PRINT GAMMA GIN BETA [FORMAT=10.8] TAUSQ; 
end;