data new; set fets.constr; 
 pi = constant("pi"); * Creates the math constant pi = 3.1415926.... ; 
 s1 = sin(2*pi*_n_/12); c1=cos(2*pi*_n_/12); * Fundamental Frequency;
 s2 = sin(2*(2*pi*_n_/12)); c2=cos(2*(2*pi*_n_/12)); * First Harmonic; 
 s3 = sin(3*(2*pi*_n_/12)); c3=cos(3*(2*pi*_n_/12)); * Second Harmonic; 
* Harmonics: Change shape while leaving period at 12; 
* Sine & Cosine: Allow shift and amplitude adjustments to match data;
* 2*pi*t/m goes through one cycle every m observations (t=1,2,...,m); 
symbol1 v=none i=join c=red; 
symbol2 v=none i=join c=black; 
Proc Gplot;
  plot contracts*date;
  plot (s1 c1)*date/overlay; 
  plot (s2 s1)*date/overlay; 

Proc Reg; 
  model contrcts = date  s1 c1;  * try adding harmonics; 
  output out=outreg predicted=p; 
proc gplot; plot (contrcts p)*date/overlay legend; 
run;