*
* SERIESF.DAT example from pp 247-255
*
open data seriesf.dat
calendar(m) 1995
data(format=free,org=columns) 1995:1 1998:12 demand
*
* Compute the trend rate, using STATISTICS to get the mean of the year over year
* difference of demand.
*
set d12 = demand-demand{12}
stats d12
compute at=%mean/12
*
* Subtract out the trend, and take the mean out of of what's left to get the
* detrended data.
*
set detrend = demand-at*t
diff(center) detrend
compute a0=%mean
disp "Trend Line is" a0 "+" at "t"
set trend = a0+at*t
*
* Do two term Fourier regression on the detrended data
*
set sin12 = sin(2*%pi/12.0*t)
set cos12 = cos(2*%pi/12.0*t)
linreg detrend
# cos12 sin12
prj seasonal
set fitted = trend+seasonal
*
* Graphs in figure 6A-6
*
graph(footer="Actual demand and trend component") 2
# demand
# trend
graph(footer="Two-term Seasonal component")
# seasonal
graph(footer="Figure 6A-6 Four-term FSA model") 3
# demand
# trend
# fitted
*
* Full seasonal model
*
set sin12   = sin(2*%pi/12.0*t)
set cos12   = cos(2*%pi/12.0*t)
set sin6    = sin(4*%pi/12.0*t)
set cos6    = cos(4*%pi/12.0*t)
set sin4    = sin(6*%pi/12.0*t)
set cos4    = cos(6*%pi/12.0*t)
set sin3    = sin(8*%pi/12.0*t)
set cos3    = cos(8*%pi/12.0*t)
set sin12_5 = sin(10*%pi/12.0*t)
set cos12_5 = cos(10*%pi/12.0*t)
set cos2    = cos(12*%pi/12.0*t)
*
linreg detrend
# cos12 sin12 cos6 sin6 cos4 sin4 cos3 sin3 cos12_5 sin12_5 cos2
*
* This computes the sequence of regressions and displays the amplitude and
* amplitude relative to the RSE. It actually isn't necessary to run the separate
* regressions to do this because the regressors are "orthogonal", that is, the
* regression X'X matrix is zero off the diagonal. When this is the case, adding
* regressors has no effect on the estimates of the earlier regressors. However,
* it's easier to just do the regressions since LINREG does most of the
* bookkeeping.
*
do i=1,5
   linreg(noprint,entries=2*i) detrend
   # cos12 sin12 cos6 sin6 cos4 sin4 cos3 sin3 cos12_5 sin12_5 cos2
   compute a=sqrt(%beta(2*i-1)**2+%beta(2*i)**2)
   disp i a a/sqrt(%seesq)
end do i
*
* Redo the regression just with the first four terms
*
linreg detrend
# cos12 sin12 cos6 sin6
prj seasonal
set fitted = trend + seasonal
set error  = demand-fitted
@uforeerrors demand fitted
*
print(picture="*.###") / demand trend seasonal fitted error