***  First run River.sas program ****;

data neuse; set neuse;
t+1;
group=floor(4*(t-1)/365)+1;    ** 4 groups per year (quarters) **;
season = 1 + mod((group-1),4);
goptions reset=all;
proc gplot data=neuse;
 plot lgold*date lkins*date/overlay legend;
symbol1 v=none i=join c=red;
symbol2 v=none i=join c=blue;

proc gplot data=neuse;
 plot lgold*date=season;
 symbol1 v=plus h=.25 i=none c=red;
 symbol2 v=plus h=.25 i=none c=cyan;
 symbol3 v=plus h=.25 i=none c=black;
 symbol4 v=plus h=.25 i=none c=green;
run;


proc arima data=neuse;
i var =lgold(1) outcov=gold;
  by group season;
  proc arima data=neuse;
  i var=lgold(1);
  proc print;
run;


data gold; set gold; nlag=-1*lag;
  lags = lag+.1*(season-1);             *** separate seasons ***;
   seas=season+.0001*ranuni(1827655);     * jitter **;
   lb=.; ub=.;
   if group=1 then do; lb=-1.96*stderr; ub=1.96*stderr;
    lbmn=-1.96*stderr/10; ubmn=1.96*stderr/10;
                                         ** 25 years = 100 quarters **;
   end;
 color="green";
 if season=1 then color="red";
 if season=2 then color="blue";
 if season=3 then color="black";

proc g3d data=gold;
scatter group*nlag=corr/
  rotate=80 tilt=80 color=color shape="cylinder" size=.4;
plot group*nlag=corr/
  rotate=80 tilt=80;
plot lag*group=corr/
  rotate=80 tilt=80;
scatter  seas*nlag=corr/
  rotate=80 tilt=80 color=color shape="balloon" size=.75;
where nlag>-4;
run;

proc gplot data=gold; plot corr*lags=season;
 symbol1 v=dot i=none c=red;
 symbol2 v=dot i=none c=cyan;
 symbol3 v=dot i=none c=black;
 symbol4 v=dot i=none c=green;
where lags < 10;
run;

 proc gplot data=gold; plot corr*lag=2 (ub lb)*lag = 1/overlay vref=0;
 symbol1 v=none i=join c=red width=3;
 symbol2 v=none i=boxt c=black;

 proc gplot data=gold; plot  invcorr*lag=2/overlay vref=-.21 -.021 0 .021 .21;
 proc gplot data=gold; plot partcorr*lag=2/overlay vref=-.21 -.021 0 .021 .21;
                     ** .21 = 1/sqrt(92 days), .021=.21/sqrt(100 quarters)**;

data gold1; set gold; if lag < 5;
proc sort data=gold1; by lag;
proc glm data=gold1;
  class season lag;
    model corr =  season/solution;
    means season/tukey;
 by lag;
run;


************ crosscorrelations ****************;

proc arima data=neuse;
    i var=lgold(1) noprint; e p=1 q=1 ml noprint;
    i var=lkins(1) crosscor=lgold(1) noprint outcov=cross;
    by group;
data cross; set cross; where crossvar="LGOLD";
  season = 1 + mod(group-1,4);
  lags = lag + .1*(season-2.5);
proc print; where corr>.8;

proc gplot data=cross;
 plot corr*lags=season;
 symbol1 v=dot i=none c=red;
 symbol2 v=dot i=none c=cyan;
 symbol3 v=dot i=none c=black;
 symbol4 v=dot i=none c=green;
    where abs(lag) < 8;
proc gplot data=cross; plot corr*lags=season/ vref=0;
 symbol1 v=none i=boxt c=red;
 symbol2 v=none i=boxt c=cyan;
 symbol3 v=none i=boxt c=black;
 symbol4 v=none i=boxt c=green;
proc gplot data=cross; plot corr*lags=season/ vref=0;
    where abs(lag) < 8;
run;