%let n = 15000; 

data future; format date date7.;
mu=7.3838;  s1b=0.9528; c1b=0.2498; a=-1.3267; b=0.0646; d=1.5278;  
sig= sqrt(6.1812/392);  sig=0.03;
date = '04nov71'd;  omega=8*atan(1)/365; 
lag2 = 9.02642 - (mu + s1b*sin(omega*(date-1)) + c1b*cos(omega*(date-1))); 
lag1 = 9.02160 - (mu + s1b*sin(omega*date) + c1b*cos(omega*date)); 
array K(&n); array LK(&n); array LLK(&n); 
do t=1 to 200; 
date = intnx('day','04nov71'd,t);
do ex = 1 to &n; 
if t=1 then do; LK(ex)=lag1; LLK(ex)=lag2; end;
L = a + B*LLK(ex); rho = (exp(L)-1)/(exp(L)+1); 
if rho> -.5279 then rho=-.5279;
trend =  mu + s1b*sin(omega*date) + c1b*cos(omega*date);
r =  D*LK(ex) + rho*LLK(ex) + sig*normal(1234567); 
K(ex)=trend+R; 
LLK(ex)=LK(ex); LK(ex)=r; 
end; output; end; 
* proc print; run; 
proc gplot; plot (trend K1-K10)*date/overlay; 
symbol1 v=none i=join c=black r=10 l=1 w=1; 
proc gplot; 
plot trend*date; 

proc transpose data=future; var K1-K&n; 
proc means noprint; 
var col1-col200; 
output out=out1 P5=P5X1-P5X200 mean=mean1-mean200 P95=P95X1-P95X200;
proc transpose data=out1 out=p5; var P5X1-P5X200; 
proc transpose data=out1 out=p95(rename=(col1=P95)); var P95X1-P95X200; 
proc transpose data=out1 out=mn(rename=(col1=mean)); var mean1-mean200; 

****** Run NEUSE1.sas first to get dataset LAST *******; 
data nonlin(keep=date forenl); set last; if date > '04nov71'd;
data id; set future(keep=date trend); 
data all; merge P5 mn P95 nonlin id; t+1;
proc print; run;  
proc print data=last( firstobs=399 obs=405); run; 
axis1 label=(angle=90 font=centb H=1.6); 
proc gplot; plot(col1 mean P95 forenl trend)*date/overlay vaxis=axis1; 
label col1= "Forecasts"; where t<101; 
title h=1.5 "Nonlinear ln(Flow) Forecasts"; 
symbol1 v=none i=join c=cyan r=3 l=3 w=3; 
symbol2 v=none i=join c=red w=3; 
symbol3 v=none i=join c=black w=3; 
run; 

data next; set river all; 
mu=7.3838; s1b=0.9528; c1b=0.2498;omega=8*atan(1)/365;
trend =  mu + s1b*sin(omega*date) + c1b*cos(omega*date);
proc gplot; plot(col1 mean P95 forenl trend lkins)*date/overlay vaxis=axis1;
label col1 = "Flow"; 
run;