MODELING PHOSPHORUS DYNAMICS
IN SHAGAWA LAKE
John VanSickle and David P. Larsen
CERL - 049
October 1978
Corval1is Environmental Research Laboratory
U.S. Environmental Protection Agency
200 S.W. 35th Street
Corvallis, Oregon 97330
PDF file saved at:
G\HOOD\LIBRARY\Digital EPA and others\EPA 0440-OCLC 18506856
10/1/2009
1

-------
epptw1®
I. Introduction
Most of the modeling done for ecological systems falls into one of two
categories, each of which has a serious disadvantage. Large systems models
contain numerous variables and parameters in attempts to maximize mechanistic
realism. Because of their complexity and large parameter spaces, these
models are difficult to corroborate with typically small ecological data
bases. On the other hand, purely statistical models summarize data well but
have little predictive value when ecosystem inputs or internal structures
change. Occasionally, however, one studies a system whose dynamics can be
described by a fairly simple mechanistic model and whose data base is large
enough to permit statistical evaluations of model adequacy.
We encountered this situation in attempting to model the dynamics of
phosphorus in a eutrophic lake. This paper describes the development and
parameter estimation of a mass balance model for phosphorus. Sufficient data
exists to statistically test the adequacy of the model, and we compare model
projections with further data when phosphorus inputs to the lake are signifi-
cantly reduced.
2

-------
II. Background
Shagawa Lake, in northeastern Minnesota, has been intensively studied by
the U.S. Environmental Protection Agency for many years. The lake has a long
history of cultural eutrophication; in 1973, a tertiary wastewater treatment
plant, funded by EPA and the City of Ely, MN, began operation in an attempt
to reduce the supply of the critical nutrients that were considered responsible
for recurrent algal blooms. Specifically, the plant was designed to remove a
large proportion of the phosphorus (P) from wastewater effluent which fed
into the lake.
Since 1973, the treatment plant has successfully kept supplies of P to
Shagawa Lake at ~20% of their former levels. It was expected that, with
reduced inputs of P to the lake, P concentrations would gradually decline as
P was lost to outflow and to deposition into lake sediments, thus relieving
the eutrophic condition.
A central objective of the continuing Shagawa study is to predict the
dynamics of P levels in the lake. We would like to predict whether or not a
steady-state condition will eventually occur if the reduction in P inputs is
maintained. Assuming a steady state will be reached, we would also like to
predict the length of time required to reach steady state, and whether steady-
state P levels will be low enough to significantly reduce algal growth.
Detailed accounts of the wastewater treatment project and the limnology
of the lake can be found in the following: Malueg et al. (1975); Larsen and
Malueg (1976); Schults et al. (1976). Summaries of the dynamics of phosphorus
and chlorophyll a in Shagawa are presented in Larsen et al. (1975) and Larsen
et al. (1978a).
ubraut
Corv*(lu R	A
200 S \Y '

-------
In conjunction with the long-term monitoring of chemical and biological
characteristics of the lake, EPA researchers have worked with several models
to describe Shagawa's history and predict its future trophic status. These
models and their predictions have been fully described by Larsen and Mercier
(1976).
The modeling effort described in this paper was, to some extent, motivated
by Hsu's recent model study of Shagawa Lake (1976). Hsu treated the phosphorus
(P) levels in the lake, recorded weekly between 1971 and 1975, as a time
series. He recognized the recurring seasonal fluctuations seen in the pre-
treatment data (Fig. 1) and described this pattern with two models. The
first model was a stochastic auto-regressive moving-average equation (Box and
Jenkins, 1970) which produces the recurrent seasonal pattern by requiring the
P level at any given week in a year to be highly correlated with the P level
in the same week of the previous year. In the second model, Hsu fit a deter-
ministic model, consisting of piecewise linear and sinusoidal functions, to
the pretreatment data, and then described the residual differences between
model and data with a Box-Jenkins moving average model.
Both models gave an excellent fit to the pretreatment pattern of P
concentrations. However, because they were both strictly empirical and based
on 1971-1972 data, neither model could predict the downward trend in lake P
levels which began following the start of wastewater treatment in 1973. Hsu
was able to follow this trend by superimposing another deterministic component
on one of the previous models, but the basic model inadequacy still remained.
Purely empirical models cannot predict novel structural changes in a time
series.
4

-------
We recognized this problem in the empirical approach and decided to try
a simple mechanistic deterministic model to describe the underlying seasonal
fluctuations in lake P levels. However, we also saw the value of the time
series techniques advocated by Hsu to deal with the residual differences
between a deterministic model and data. We used time series methods to study
the adequacy of a mechanistic P model and to provide confidence limits for
the P trajectories which it would predict.
This paper describes the development of a mechanistic model and treats
in some detail the estimation of model parameters, analysis of residuals, and
generation of confidence intervals. The actual model predictions as compared
with observed lake P levels and the resulting implications for recovery of
the lake are discussed in another paper (Larsen et al., 1978a). Here the
emphasis will be on the analytical methods which generated the model predic-
tions and the limitations on these methods imposed by the data.
III. The Deterministic Model
We used a mass balance equation to describe changes in lake P. In the
model, P is dealt with as a concentration by dividing mass flow rates by the
lake volume. As illustrated by Figure 2, the change in P levels is a sum of
inflow, outflow, deposition and sediment release rates.
- v - [P] ¦ °pcp] + v	0)
where: [P] = Concentration of total P, (pfl/l)
J = Rate of P supply from external sources
V = Lake volume
Q = Rate of water outflow
5

-------
Op = Instantaneous loss rate of P to the lake sediments.
R = Rate of P release from the sediments
With the exception of the release terra, R/v> equation (1) is an accepted
description of an instantly well-mixed lake system which loses P to the
sediments (Vollenweider, 1975; Sonzogni et al., 1976; Larsen et al., 1978a).
Sbagawa Lake is unusual in that its P dynamics have been dominated by a
significant feedback of P from the sediments. This release of P appears in
Figure 1 as the sudden increases which begin about the 10th and 25th weeks of
each year. Larsen et al. (1978b) discuss the mechanisms for this phenomenon,
and their analysis indicates that both the spring and summer releases begin
suddenly, proceed at a constant rate for a number of weeks and then end
suddenly. Thus, we modeled the release rate as
fRi, 11 ^ t ^ tg
R2, t3 < t < t4	(2)
0, otherwise
The switch times tx - t4 refer to specific weeks of an arbitrary year,
so R(t) is a periodic function having a period of 1 year. A final assumption
in equation (1) is that we take 0p to have one of two constant values depending
on the season of the year, viz.
fa s if lake is ice-free (week 17 to week 46)
v* ) p
f o w if lake is ice-covered (week 46 to week 17)
v P
This assumption is more in agreement with the mechanics of P deposition than
the use of a fixed ay
6

-------
Equation (1) with a piecewise constant op and R(t) specified by (2),
constitutes the deterministic model for P in Shagawa Lake. The schedules for
P input (J) and water flow rate (Q) are provided on a weekly basis from flow
and P loading data whose measurement has been described elsewhere (Malueg et
a]., 1975; Larsen et air, 1975; Larsen and Malueg, 1976).
IV. Preliminary Model Analysis
Because equation (1) is linear we can address the questions of model
steady-state behavior which we previously raised in section II with respect
to actual lake dynamics. Suppose the model is started at an arbitrary initial
condition. Further, suppose the same weekly schedules for J and Q are used
in the model year after year so that all parameters and forcing functions are
periodic with a period of 1 year. Then there exists a unique steady-state
for [P(t)] which has the same period as that of the forcing function, J(t)/V
+ R(t), i.e. the steady state has a period of 1 year (Yakubovich and Starzhin-
skii, 1975). This result is true as long as the unforced model, with J(t)/V
+ R(t) = 0, has no periodic steady state. By a result of Floquet theory
(D1Angelo, 1970), the unforced model can have a periodic steady state only if
state is periodic in the forced model.
This steady-state result was borne out by simulations in which the
initial condition for the model was taken to be the value of [P] at the end
of 1972. Following parameter estimation, the model was run with a reduced
+ a dt = 0 (D'Angelo, 1970, p. 194)
P
(3)
Since Q and op are not negative, the integral cannot be zero and the steady
7

-------
schedule for Inputs typical of a post-treatment year and two different water
flow schedules - one for a "wet" year and one for a "dry" year. In both
cases, steady-state was reached within three years.
Thus, if the model is valid, we would expect to see the lake P levels
show a new steady-state type of behavior within a few years after P inputs
are reduced. Of course, an exact steady state will not be seen because of
year-to-year variations in the loading and washout rates, J and Q.
V. Parameter Estimation
The next stage in the model development is the estimation of parameters
for equation (1) based on the pretreatment data time series of Figure 1.
For the estimation, the pretreatment interval from the onset of ice cover in
late 1971 until the onset of ice cover in 1972 was partitioned into an ice-
covered and an ice-free period. The 1971-1972 ice-covered interval (23
weeks) was used to estimate cjpw, Rlt tj and t2; the ice-free period during
1972 (29 weeks) was used to estimate a^s, R2» t3 and t4. The remainder of
the pretreatment period, made up of the first 46 weeks of 1971 and the last 6
weeks of 1972, was used to corroborate the model (Fig. 3).
We used a least-square criterion for the estimation, i.e.
Minimize: C = 2	^
with respect to the parameters OpW, a^s, Ri R2» tj. t2 t3 t4. In equation
(3). [P]|, and [P]. are weekly concentrations from the data and model, respec-
K	IN
tively, and the week index k runs over the appropriate estimation interval.
The large size of the parameter space strongly suggests using a directed
search technique to minimize C. There are several gradient search algorithms
8

-------
available for this type of optimization (Pierre, 1969), and we decided to use
a gradient approach. However, the problem of minimizing C has an unusual
feature — the gradient of C must be computed with respect to the switch time
parameters, tj - t4.
Fortunately, methods which have been derived in the context of optimal
control theory for calculating this type of gradient do exist. It is not
difficult to recast the minimization problem (3) in control theoretic terms.
We view equation (3) as a cost functional for a linear tracking problem with
fixed initial and final times (Kirk, 1970).
In the remainder of this section we will sketch the results from the
calculus of variations which compute the gradient of C. We follow closely
the method of Hasdorff (1976), which should be consulted for derivations and
other mathematical details. One preliminary note: the cost functional C
will be restated as an integral to simplify the mathematics. Implementation
of the results requires that discrete time be used, i.e., C takes the form of
equation (3) and equation (1) is restated as a difference equation with an
update interval of one week. However, since equation (1) is linear with rate
coefficients having small absolute values, discretization of the continuous
equations is straightforward and details will be omitted.
Let v = [ct w, Rx t1( t2]+, where the + denotes vector transpose and
P
~ signifies a vector. For estimation of these parameters we restate C as
where [tQ, tf] is the 1971-1972 ice-covered period, and where [P] has been
renamed xx. We also restate equation (1) as
(4)
9

-------
dx,
*	Ai A,
It = fl(x' V)
with fi given by the right-hand side of equation (1) and x = [xlf x2]+.
Define the state variable x2(t) by
dx2	*
-ft = f2(x, V) = (Xj - Xi)2, x2(t0) = 0	(5)
This definition allows us to write the cost functional in terms of a system
state variable, viz.
C = 0 (x(tf)) = x2(tf)	(6)
The notation combines the systems dynamical equation and the cost functional
into a single state equation
| = ra, ?)	(?)
Finally, we need to define adjoint variables, X (t), also known as costate
variables or Lagrange multipliers (Kirk, 1970). The adjoint variables satisfy
the state equation
= -T X : X(tf) = V 0(x(t »	(8)
at	x	i a i
In equation (8), >x is the matrix of partial derivatives whose ij th element
is given by

10

-------
Also, Vx signifies a gradient with respect to the vector x, i.e.
p>«x(tf))
»0(x(tf))
8x5
Armed with these definitions, Hasdorff uses variational methods to calculate
the gradient of an arbitrary cost functional, 0(x(tf)), with respect to
several forms of initial conditions and input and control parameters. For
the system described by equations (4) and (5) the gradient is
9(apw, Rlt tlf t2) =
£ -
t *(t)dt
fit
f* £(t)dt
rt2
K <
Tcti) [^(xctj.o) - ?(x(t1),R1)j
X+(t2) pCxCta)^!) - ?(x(t2),R2)j
(9)
Each component in the gradient vector corresponds to the parameter having the
same position in the argument list of g. In addition,
K-
3fj dfs
Ba.
w'
da
w
V
afj af2
MI' br7
Calculation of the gradient is accomplished with the following sequence:
1. Given an initial guess for the parameter vector, v, the state
equation (7) is integrated numerically from tQ to tf. This calcula-
tion also gives the value of C.
11

-------
2.	The final state, xft^), provides the initial condition for the
adjoint equations via (8).
3.	The adjoint system (8) is integrated backwards in time from tf to
t and the values of %(t) are stored,
o
4.	Equation (9) uses the stored values of \(t) to compute the gradient
9-
5.	The values of C and g are passed on to a gradient search algorithm
which uses them to change the parameter vector such that C decreases.
The sequence is repeated until C has been satisfactorily minimized. For step
5 of the sequence, we used the fast and efficient Davidon-Fletcher-Powel1
conjugate descent algorithm (Pierre, 1969). The algorithm is implemented by
subroutines FMFP and DFMFP of the IBM System/360 Scientific Subroutine package.
The algorithm converged quickly for this problem, but, due to the the noise
in the data, the cost functional was fairly insensitive to small changes in
some parameters. Thus, estimation runs were made in groups, with each run of
a group having different initial guesses for the parameters.
Table 1 gives the best estimates for the model parameters. The fitted
model and data are shown in Figure 3. As measured by the sum of squared
residuals, the model showed a 25% better fit during the corroboration interval
than it did during the estimation interval. Therefore we concluded that
equations (1) and (2) with the parameter estimates of Table 1 provide an
adequate deterministic model for the P dynamics seen in Figure 1.
VI. Analysis of Residuals
A. The Stochastic Model
The process of parameter estimation produced a "best fit" »del according
12

-------
to the least squares criterion. However, this type of fitting does not
necessarily eliminate systematic departures, through time, between model and
data. We can test for these trends by analyzing the model residuals using
the time series techniques of Box and Jenkins (1970).
The residual time series ek for 1971-1972 is shown in Figure 4. The
first step in the analysis of this series is calculation of the mean and
variance. Since the mean value of the series is effectively zero, we next
studied the autocorrelation structure of e^. The theoretical autocorrelation
function at time lag j is defined by
E [(e. - MJ(ek+i - Me)3
p =	£	?	S-J	S-	(10)
J	Qf 2
e
where E is the expectation operator, j and k are discrete time indices, and
Me and are the mean and variance, respectively, of the series e^.
Figure 5 shows the function pi as estimated for the e's of Figure 4. The
«J
estimated values of p. indicate whether the	e-series is covariance stationary.
If the e-series is stationary, then pj will	approach zero rapidly with
increasing j. The dashed lines in Figure 5	show the standard error of the
estimate of p.. Estimated values of an autocorrelation which is theoretically
J
zero are approximately normally distributed. Hence values of p^, such as
Pi7» which are theoretically zero but fall outside the line have about a one-
third probability of occurring (Box and Jenkins, 1970, p. 178).
Another restriction on the reliability of pj estimates is provided by N,
the number of data points available. Here N = 104, and Box and Jenkins
(1970, p. 33) claim that estimates of p^ are unreliable for j > N/4. The
deterministic model and the lake P dynamics operate on a fundamental period
13

-------
of 1 year, so we expect to see a peak in at j = 52.	This peak does appear
in Figure 5, but the above criterion requires at least 4 years of data for
one to be statistically assured of its importance. As	it is, we cannot
depend on the estimated value of p. beyond lag 26.
J
Based on Figure 5, we decided to try a stationary model for ek with a
theoretical autocorrelation function which is effectively zero beyond lag 2.
A second order moving-average model of the following form is indicated:
ek = ak " 61 ak-l " 02 ak-2	(11>
In equation (11), the a's are independent random shocks having zero mean and
constant variance o£.
a
The model (11) was fit to the series e^. using an estimation algorithm
given in chapter 7 of Box and Jenkins (1970). The results of the estimation
produced
ek = ak + 0.55 ak-1 + 0.2 aR_2, a* = 30.6	(12)
Finally, we apply diagnostic checks to assess the adequacy of the model
(12). The procedure is to use (12) and the residuals ek to generate a
series a.. If equation (12) is adequate, the ak's should be uncorrelated
l\
with a zero mean. A white noise test is performed on the autocorrelation
function of the a's (Box and Jenkins, 1970, p. 299). This test confirmed
In
the adequacy of (12).
The preceding analysis shows that the model residuals, ek, can be
described by a stationary, random time series with a zero mean and a low
oraer of autocorrelation. We interpret this result to mean that the mass
balance model (1) adequately represents a^l of the deterministic information
about P dynamics which can be extracted from the 1971-1972 data.
14

-------
B. Forecasts and Confidence Limits
With the aid of equation (12) we can place confidence limits on the
predictions of the deterministic model (1). We assume that the nature of the
random errors, e^, remains unchanged following the start of wastewater treat-
ment in 1973. Following the method of Box and Jenkins (1970 chapter 5), we
forecast values for the series e^ and superimpose the forecast on the deter-
ministic predictions. The error variance of the forecast gives confidence
limits for the prediction of ek.
Let e. . = forecast of the e-series for j weeks ahead of the time
K+J
origin k. If we take k as the last week in 1972, then, except for the first
two weeks in 1973, the forecasts ek+j are identically zero, since pe = 0 and
the model (12) is stationary. Therefore, the expected future value of the
residuals is zero, and we center the confidence limits on the trajectory of
the deterministic model.
For j >3, the confidence intervals are defined by
195% confidence interval! _ p + n	.q	H31
at lead time j J " Pt*j 1 U°-05/2 °e	(l3)
where P. . is the deterministic prediction at time k+j, and
k+j
deviate exceeded by the fraction e/2 of the unit Normal distribution.
VII. Model Predictions
Figures 6 and 7 show model predictions compared with data for 1973 74
and 1975-76. The two predicted trajectories employ observed flow rates
during 1973-1976 for the parameters Q(t). The prediction which assumes
15

-------
treatment uses the observed P inputs for J(t), which are about 80% lower than
pretreatment inputs (Larsen et. al., 1978a). We made the prediction which
assumes no treatment by employing the average P input rate observed during
the two pretreatment years 1971-1972. Confidence intervals are shown centered
on the trajectory which assumes treatment, but intervals of equal width would
also apply to the "no treatment" curve.
Since observed values for J and Q were used, Figures 6 and 7 are not
predictions in the truest sense. To make such predictions, we would also
need to forecast time series for J and Q, both of which are functions of
river discharge, and thus of weather patterns. However, we compared simula-
tions of the model which employed Q(t) for a typical "high" water year as
opposed to a "low" water year and found that the P trajectories were never
separated by more than 4 pg/1. Based on this trial, we feel that the model
would provide good predictions for future normal water years.
The observed lake P response to reduced loading appears to be delayed
since it follows the "no treatment" prediction more closely than the "treat-
ment" prediction during most of 1973. However, by 1975 and certainly 1976,
the observed lake P levels track the "treatment" prediction quite closely and
consistently fall within the confidence intervals. Following a sharp reduction
in P inputs, equations (1), (2) and (12) clearly provided an accurate descrip-
tion of P dynamics in Shagawa once the lake response has become stabilized.
Larsen et al. (1978a) give a statistical summary of the results shown in
Figures 6 and 7 and they discuss the Figures' implications for improvement in
the trophic status of the lake.
16

-------
VIII.Summary and Conclusions
We have outlined the methodology used to develop a predictive model of
phosphorus dynamics in Shagawa Lake. The modeling process involved a sequence
of steps, beginning with selection of a simple mass balance equation and
mathematical formulation of the various rate functions in the equation. In
the next stage, we estimated model parameters using a least squares criterion
applied to pretreatment data of 1971-1972. The fitted model was then subjected
to a time series analysis of its residuals to uncover systematic departures
from the data. Upon finding no trends in the residuals, we judged the deter-
ministic model to be adequate and went on to generate confidence limits from
the stochastic model of the residuals. Finally, we made predictions of P
levels for 1973-1976 and compared them with data.
In constructing this model we have taken an approach more akin to
conventional statistics than to currently popular practices of systems
ecology. That is, we began with the constraints imposed by the data base,
constraints of sampling frequency, sampling duration, and types of variables
sampled. Our goal was to build a predictive model which could extract the
maximum amount of information about lake P dynamics from the available data,
and would require as few assumptions about mechanisms or parameter values as
possible. The recent work of Walker (1977) is an excellent discussion of
this approach applied to a variety of lake water quality models.
The data base for Shagawa Lake constitutes one of the longest, most
comprehensive records available for variables pertinent to algal-nutrient
dynamics. Thus, even though equation (1) appears fairly elementary, it may
be the most sophisticated model of its type which can be supported by currently
available data.
17

-------
References
Box, G.E.P., and G.M. Jenkins. 1970. Time Series Analysis: Forecasting and
Control. Hoi den-Day, San Francisco, 575 pp.
D'Angelo, H. 1970. Linear Time-Varying Systems: Analysis and Synthesis.
Allyn and Bacon, Boston. 346 pp.
Hasdorff, L. 1976. Gradient Optimization and Nonlinear Control. John Wiley
and Sons, New York. 264 p.
Hsu, D.A. 1976. A sutdy on phosphorus concentrations in Shagawa Lake,
Minnesota. Unpubl. manuscript. 18 p.
Kirk, D.E. 1970. Optimal Control Theory. Prentice-Hall, Englewood Cliffs.
452 pp.
Larsen, D.P., K.W. Malueg, D.W. Schults, and R.M. Brice. 1975. Response of
eutrophic Shagawa Lake, Minnesota, U.S.A., to point source phosphorus
reduction. Verh. Internat. Verein. Limnol. 19:884-892.
Larsen DP and H T Mercier. 1976. Shagawa lake recovery characteristics
as depicted by predictive modeling. U.S. Environmental Protection
Agency Report EPA-600/2-76-079. Corvallis Environmental Research Labora-
tory.
Larsen, D.P., and K.W. Malueg. 1976. Limnology of Shagawa Lake Minnesota,
prior to reduction of phosphorus loading. Hydrobiologia 50:177-189.
Larsen, D.P., J. Van Sickle, and K.W. Malueg. 1978a. The effect of wastewater
phosphorus removal on Shagawa Lake: phosphorus supplies, lake phosphorus,
and chlorophyll a. Manuscript in preparation.
Larsen DP D W Schults, and K.W. Malueg. 1978b. Summer internal phospho-
rs s^Shagawa Lake, Minnesota. Manuscript in preparation.
Maluea K W DP Larsen D.W. Schults, and H.T. Mercier. 1975. A six-year
waterf phosphorus! and nitrogen budget for Shagawa Lake, Minnesota.
Journal of Environmental Quality 4.236 til.
Pierre, D.A. 1969. Optimization Theory with Applications. Wiley, New York.
612 p.
Schults n w K W Maluea and P.O. Smith. 1976. Limnological comparison
of'cuitiiraliy'eutrophlc Shagawa Lake and adjacent oligotrophy Burts,de
Lake, Minnesota. The American Midland Natural,st 96:160-178.
Sonzogni, W.C., P.O. Uttormark, and G.F. Lee 1976 The phosphorus residence
ti«4 model: theory "><* application. Water Res. 10.429
18

-------
Vollenweider, R.A. 1975. Input-Output models with special reference to the
phosphorus loading concept in limnology. Schweiz. Z. Hydro!. 37:53-83.
Walker, W.W., Jr. 1977. Some Analytical Methods Applied to Lake Water
Quality Problems, ch. 2. Ph.D. Thesis, School of Engineering, Harvard
University. 210 pp.
0
Yakubovich. V.A. and V.M. Starzhinskii. 1975. Linear Differential Equations
with Periodic Coefficients, V.l. John Wiley & Sons, New York. 386 pp.
19

-------
TABLE 1. Parameter estimates for equations (1) and (2)
Ice-covered season	Ice-free season
II
*0*
0.035 wk-1
°p
= 0.072 Wk-1
Ri =
280 kg/wk
r2
= 500 kg/wk
=
week 12
t3
= week 25
ii
week 14
t4
= week 34
20

-------
100
TIME (weeks)
TREATMENT
BEGINS
Figure 1.	Weekly average total phosphorus levels in Shagawa Lake,
1971-1972.

-------
Figure 2.	Sources arid losses of phosphorus in Shagawa Lake.

-------
100
1971A | +1972
40
20
CORROBORATION-
I	l_

— ESTIMATION
1 1
0
10
20
30
40
50
10
20
30
40
50
TIME (weeks)
TREATMENT
BEGINS
Figure 3.
Fitted model trajectory, 1971-1972.

-------
CP
o>
LJ
O
UJ
CC
UJ
II.
U-
Q
-J
<
9
c/>
UJ
o:
-20
TIME (months)
Figure 4.
Weekly residual differences, e^-
1971-1972.
Data - [P]?0del).
24

-------
0.96 h
o
"o
£
Ts>
ui
0.72
0.48
u)
(X.
cc.
o
o
o
0.24h
-0.241—
LAG NUMBER
Standard
Error
Figure 5.	Autocorrelation of residual series.
25

-------
100
1973< | +I974
(95% CONFIDENCE
BANDS SHOWN)
TREATMENT
BEGINS
TIME (weeks)
Figure 6.
Model predictions compared with data, 1973-1974.

-------
rv>
CJ»
a.
CO
3
-------