options ls=76; 
***********************************
* Create transfer function data.  *
***********************************;
Data a; pi = 4*atan(1); 
 do t = 1 to 107; Ct = cos(2*pi*t/15); St=sin(2*pi*t/15); 
   *X = 3*Ct + 4*St; 
    X = 5*ranuni(123);
   Y = 5 + X + 2*X1 + 3*X2; 
output; X2=X1; X1=X; end; 
proc print data=a (obs=5); title "Bivariate Data";
proc gplot; plot (X Y)*t/overlay; 
symbol1 v=none i=spline c=black r=2;  label X="X Y";
title "Wave run through a Transfer Function"; 

*******************************************
* Compute Ssq by regression, PROC SPECTRA *
*******************************************;
proc reg; model Y = St Ct; 
proc spectra P adjmean out=out1; var Y; 
proc gplot data=out1; plot P_01*freq; 
 symbol1 v=diamond i=needle; 
proc spectra data=a adjmean P S cross Ph A K out=out2; 
  var X Y; 
  data out2; set out2;
  ********************************************* 
  *  Add multiples of pi to smooth everything *
  *********************************************; 
  if PH_01_02 < 0 and FREQ < 1.5 then PH_01_02 = PH_01_02 + 4*atan(1) ; 
  if PH_01_02 < 0 then PH_01_02 = PH_01_02 + 8*atan(1); 
  if round(freq,0.0001) = 2.5731 then PH_01_02 = PH_01_02 + 4*atan(1); 
proc print; run;
proc gplot data=out2; plot PH_01_02*freq; 
symbol1 v=none i=join; run;