options nodate nonumber ps=150 ls=80;

title 'RPD -- Sec 8.2 -- Quarterly Beer Production';

data beer;
input time production;
datalines;
1      36.14
2      44.60
3      44.15
4      35.72
5      36.19
6      44.63
7      46.95
8      36.90
9      39.66
10     49.72
11     44.49
12     36.54
13     41.44
14     49.07
15     48.98
16     39.59
17     44.29
18     50.09
19     48.42
20     41.39
21     46.11
22     53.44
23     53.00
24     42.52
25     44.61
26     55.18
27     52.24
28     41.66
29     47.84
30     54.27
31     52.31
32     42.03
;


symbol interpol=join   /* connect the points */
       value=dot;

proc gplot;
plot production*time;
run;


proc iml;

use beer;
read all;

Y = production;

   /* Pure annual-cycle second-degree trigonometric model (ignoring trend): */

X1 = j(nrow(Y),1)||cos(2*constant('PI')*time/4)||
                   sin(2*constant('PI')*time/4)||
                   cos(2*constant('PI')*2*time/4)||
                   sin(2*constant('PI')*2*time/4);

reset fuzz;

print X1;

X1 = X1[,1:4];   /* use only first four columns --- last column is all zeros */

XTX1 = X1`*X1;

XTXINV1 = inv(XTX1);

BETAHAT1 = XTXINV1*(X1`*Y);

SSRES1 = Y`*Y - BETAHAT1`*X1`*Y;

S2_1 = SSRES1/(nrow(X1)-ncol(X1));

print X1 XTX1 BETAHAT1 S2_1;

    
   /* Model augmented with a linear trend: */

X2 = X1||time;

XTX2 = X2`*X2;

XTXINV2 = inv(XTX2);

BETAHAT2 = XTXINV2*(X2`*Y);

SSRES2 = Y`*Y - BETAHAT2`*X2`*Y;

S2_2 = SSRES2/(nrow(X2)-ncol(X2));

print X2 XTX2 BETAHAT2 S2_2;


   /* To test whether first-degree model is adequate ... */

T0 = BETAHAT2[4]/sqrt(XTXINV2[4,4]*S2_2);

PVALT0 = 2*(1-probt(abs(T0),nrow(X2)-ncol(X2)));

print T0 PVALT0;


   /* To test if any trigonometric terms are needed at all ... */

X_R = j(nrow(Y),1)||time;

XTXINV_R = inv(X_R`*X_R);

BETAHAT_R = XTXINV_R*(X_R`*Y);

SSRES_R = Y`*Y - BETAHAT_R`*X_R`*Y;

F0 = ((SSRES_R-SSRES2)/3)/S2_2;

PVALF0 = 1-probf(F0,3,nrow(X2)-ncol(X2));

print F0 PVALF0;


   /* Compute predicted values for final model: */

X = j(nrow(Y),1)||cos(2*constant('PI')*time/4)||
                  sin(2*constant('PI')*time/4)||
                  time;

XTXINV = inv(X`*X);

BETAHAT = XTXINV*(X`*Y);

YHAT = X*BETAHAT;

   
   /* Save predicted values to a SAS data set: */

create fitted from YHAT [colname='yhat'];
append from YHAT;

quit;


data beer;
merge beer fitted;   /* combine data sets 'beer' and 'fitted' */
run;


proc gplot;
plot (production yhat)*time / overlay legend;   
   /* overlay plots of data and fitted values, and provide a legend */
run;