%let n=48;  *n rows -> n+2 observations ; 
%let r = 2000; * reps of each test; 
options symbolgen;

proc iml; 
  Y = shape(0, &n, 3); 
  X = shape(0,&n,2); 
  FIRST = {0,0}//Y[ ,1]; 
  df0=&n; df=df0-1; df2=df0-2; 
  do t = 1 to &n; 
   X[t,1]=1; X[t,2]= t - (&n+1)/2; end; 
  XP2 = X*inv(X`*X)*X`; 
  XP1 = X[ ,1]*inv( X[ ,1]`*X[ ,1] ) * X[ , 1]`; 
  outdat = shape(0,1,10); 
  roots = {0, .3, .6, .8, 1 }; 
  zoots = {.7, .8, .9, .95, .98,  1}; 
  labs={"root1" "root2" "rho" "tau" "rho_M" "tau_M" "F1" "rho_t" "tau_t" "F2"}; 

  **** Simulation *****; 
  do i = 1 to 5;  ** root 1; 
  do j = 1 to 6;  ** root 2;
  do rep = 1 to &r;  ** reps; 

  in=(i=1)*(j=1)*(rep = 1); ** first series out as check ; 

  a1 = (roots[i,1]+zoots[j,1]); ** AR coefficents ; 
  a2 = -1*roots[i,1]*zoots[j,1];
  ** initialize Y ; 
  Y[1,3] = normal(6545387);Y[1,2]=a1*Y[1,3] + normal(6548277); 
  Y[1,1] = a1*Y[1,2] + a2*Y[1,3] + normal(2637165);

       if in then do; first[1,1] = Y[1,3]; ** first series out; 
       first[2,1]=Y[1,2]; first[3,1]=Y[1,1]; end; 

  do t = 2 to &n; ** run 2 through n; 
    Y[t,3]=Y[t-1,2]; Y[t,2]=Y[t-1,1]; 
    Y[t,1]=a1*Y[t,2]+a2*Y[t,3]+normal(6577655);  
	if in then first[t+2,1]=Y[t,1];
  end; 

  ** Regress out deterministic effects ; 
  R2 =Y-XP2*Y; 
  R1 = Y - XP1*Y; 
  * print R1 R2; 

  ** Mean 0 tests ; 
  YY = Y[ ,1]-Y[ ,2]; 
  Z = Y[ ,2:3]; Z[ ,2]=Z[ ,1]-Z[ ,2]; 
  COVB=INV(Z`*Z);
  b = COVB*Z`*YY; 
  YPY = YY`*YY; 
  SSE = (YPY - B`*Z`*YY); 
  MSE = (1/(&n-2))*SSE; 
  COVB = COVB*MSE;  
   * RHO ; 
  OUTDAT[1,1]=ROOTS[I,1]; OUTDAT[1,2]=ZOOTS[J,1]; 
  OUTDAT[1,3] = &n*B[1,1]*INV(1-B[2,1]);
   * TAU ; 
  OUTDAT[1,4] = B[1,1]/SQRT(COVB[1,1]); 
 
** Estimated mean ;  
  YY = R1[ ,1]-R1[ ,2]; 
  Z = R1[ ,2:3]; Z[ ,2]=Z[ ,1]-Z[ ,2]; 
  COVB=INV(Z`*Z);
  b = COVB*Z`*YY; 
  YPY = YY`*YY; 
  SSERED= YPY - B`*Z`*YY ;
  MSE = (1/(&N-3))*SSERED; 
  COVB = COVB*MSE; 
   * RHOMU ; 
  OUTDAT[1,5] = &n*B[1,1]*INV(1-B[2,1]);
   * TAUMU ; 
  OUTDAT[1,6] = B[1,1]/SQRT(COVB[1,1]);
   * F1 ; 
  OUTDAT[1,7] = (SSE-SSERED)/MSE; 

** Detrended tests;  
  YY = R2[ ,1]-R2[ ,2]; 
  Z = R2[ ,2:3]; Z[ ,2]=Z[ ,1]-Z[ ,2]; 
  COVB=INV(Z`*Z);
  b = COVB*Z`*YY; 
  YPY = YY`*YY; 
  SSERED = YPY - B`*Z`*YY;
  MSE = (1/(&N-4))*SSERED; 
  COVB = COVB*MSE; 
   * RHO_t ; 
  OUTDAT[1,8] = &n*B[1,1]*INV(1-B[2,1]);
   * TAU_t ; 
  OUTDAT[1,9] = B[1,1]/SQRT(COVB[1,1]);
   * F2 ; 
  OUTDAT[1,10] = (SSE-SSERED)/(2*MSE); 

** Append results and send to dataset; 
if in then send=outdat; else send = send // outdat; 
END; END; END; 
* print send [COLNAME=LABS];  
* print first;
CREATE SIMDAT FROM SEND [COLNAME=LABS]; 
APPEND FROM SEND; 

PROC PRINT DATA=SIMDAT(obs=200); RUN;  
data simdatall; set simdat; array rej(8);   
rej1 = (tau<-1.95); rej2 = (rho<-8.1); rej3=(tau_M<-2.86);
rej4= (rho_M<-14.1); rej6 = (tau_t<-3.41); rej7=(rho_t<-21.7);

/*  Code to check the calculations on first replicate ; 
create dataout from first; 
append from first; 


data in; set dataout (rename=(col1=Y)); t+1; 
proc arima data=in; i var=Y stationarity=(adf=(1)); 
data next; set in; D = Y-Y1; d1=Y1-Y2; 
output; y2=y1;y1=y; retain; 
proc reg; 
  model D = D1 t Y1/SS1;
    test t=0, Y1=0; 
    test intercept=0, t=0, y1=0;
  model D = D1 Y1/ss1; 
    test Y1=0; 
    test intercept=0, Y1=0;  
  model D = Y1 D1/noint; 
  * For RHO, SAS seems to use n=number of obs in regression, divides lag level
  coefficent by (1-augmenting difference coefficient) ; 
*  */
run; 

proc means data=simdatall noprint; var tau tau_m tau_t rej1-rej4 rej6 rej7; 
class root1 root2; 
output out=outpow mean = mn0 mn1 mn2 ptau prho ptau_m prho_M ptau_t prho_t;
proc means data=outpow; var ptau--prho_t;
proc print data=outpow; 
data outpow2; set outpow; keep root1 root2 Y0 Y_M Y_t type cval; 
Y0 = prho; Y_M=prho_M; Y_t=prho_t; type="Rho"; cval = "blue"; output; 
Y0 = ptau; Y_M=ptau_M; Y_t=ptau_t; type="Tau"; 
root1=root1 - .03; root2=root2+.01;
cval="red"; output; 
proc g3d data=outpow2; 
scatter root1*root2=Y0/ zmax=1 color=cval shape="prism"; 
scatter root1*root2=Y_M/ zmax=1 color=cval shape="prism"; 
scatter root1*root2=Y_T/ zmax=1 color=cval shape="prism"; 
title h=1.5 c=blue f=centB "Rho " c=red " Tau"; 
label Y0 = " 0"; label Y_M = "Mu"; label Y_t ="t";
run;