# Removing observational noise from fisheries-independent time series data using ARIMA models.

Abstract--Abundance indices derived from fishery-independent
surveys typically exhibit much higher interannual variability than is
consistent with the within-survey variance or the life history of a
species. This extra variability is essentially observation noise (i.e.
measurement error); it probably reflects environmentally driven factors
that affect catchability over time. Unfortunately, high observation
noise reduces the ability to detect important changes in the underlying
population abundance. In our study, a noise-reduction technique for
uncorrelated observation noise that is based on autoregressive
integrated moving average (ARIMA) time series modeling is investigated.
The approach is applied to 18 time series of finfish abundance, which
were derived from trawl survey data from the U.S. northeast continental
shelf. Although the a priori assumption of a
random-walk-plus-uncorrelated-noise model generally yielded a smoothed
result that is pleasing to the eye, we recommend that the most
appropriate ARIMA model be identified for the observed time series if
the smoothed time series will be used for further analysis of the
population dynamics of a species.

**********

Time series of species abundance from fishery-independent surveys, such as bottom trawl or acoustic surveys, are important in monitoring temporal change in the abundance of marine populations. For commercially important species, catch and effort data from the commercial fishery may be available, allowing estimation of temporal trends of the stock population by means of stock assessment models (e.g., virtual population analysis). However, such records are not available for many species, especially those with little commercial (but perhaps significant ecological) value. Fishery-independent surveys may thus constitute the only source of information for assessing temporal changes in the abundance of these species (Pennington, 1985; Helser and Hayes, 1995).

Annual estimates of abundance derived from fisheries-independent surveys are typically regarded as providing a relative measure of population abundance (i.e., they are indices of abundance, not true estimates of total population size) (Grosslein, 1969; Clark, 1979). Thus, the expected value of the abundance index (e.g., mean catch-per-tow for trawl surveys) is regarded as proportional to the size of the actual population, although the constant of proportionality (the catchability) is unknown. As such, relative changes in an abundance index should reflect similar relative changes in the actual population, and trends in the time series of such an index should reflect similar trends in the corresponding population.

Unfortunately, abundance indices derived from large-scale fishery-independent surveys typically exhibit interannual variability much higher than one would expect from withinsurvey variance (Byrne et al., 1981; Pennington, 1985). Part of the variability in such indices is presumably due to the variability in the underlying population--a variability that is caused by population-dynamic processes such as recruitment. However, part of the variability is due to observation noise that arises from both within-survey sampling variability because of the heterogeneous distribution of many fish stocks (Byrne et al., 1981), and because of environmentally driven factors that affect catchability over time (Byrne et al., 1981; Collie and Sissenwine, 1983). Low signal-to-noise ratios in abundance indices that are due to high observation noise reduce chances of detecting important changes or trends in actual population abundance. Variability due to within-survey sampling can be reduced (before the fact) by adding more stations to a survey, but additional stations will not reduce variability caused by changes in catchability.

Time series modeling using autoregressive integrated moving average (ARIMA; Box and Jenkins, 1976) models provides an approach to removing observation noise from abundance estimates. ARIMA models are frequently used in economic forecasting (Enders, 2004) and are becoming more common in fisheries research. Recent applications of ARIMA models to other fisheries problems include forecasting monthly landings in the Mediterranean (Lloret et al., 2000), testing theories of population dynamics (Becerra-Munoz et al., 1999), and modeling nutrient dynamics in an upwelling system (Nogueira et al., 1998).

In the context of reducing the influence of observational noise in time series data, Cleveland and Tiao (1976) first developed a noise-reduction and smoothing algorithm for processes that could be described by an ARIMA time series model. Their approach requires that the ARIMA model for the unobserved, underlying process be known. This known model, in turn, uniquely determines the ARIMA model for the observed time series contaminated by observation noise and allows one to estimate the variance of the observation noise. Unfortunately, although the ARIMA model for an unobserved, underlying process may be known in some instances (from theory, perhaps), in many cases the model for the unobserved process will be unknown.

Box et al. (1978) extended Cleveland and Tiao's (1976) ideas and developed a noise-reduction algorithm based on the ARIMA model for the observed time series. However, the ARIMA model for the observed time series merely constrains, but does not determine, the model for the unobserved, underlying process; it provides only an upper bound for the observation error variance. Consequently, this approach generally requires an external estimate of the observation error variance to determine the appropriate level of noise reduction.

Pennington (1985) first applied these ARIMA-based time series modeling techniques to smoothing abundance indices derived from trawl survey data. He assumed that an observed abundance time series reflected a combination of the underlying population abundance and independent, uncorrelated, and multiplicative observation noise (the latter arising perhaps from environmentally driven changes in catchability). He further assumed that both the (log-transformed) observed time series and unobserved population process could be represented by ARIMA models. Pennington (1985) then developed an alternative algorithm to that of Box et al. (1978); his derivation allowed particular simplification in the case where the underlying population process could be modeled as a random walk. In this simple case, the resulting noise reduction filter is an exponentially-weighted average of the observed time series for the endpoint of the time series (Pennington, 1985). More importantly, the observation error variance can be easily estimated from the ARIMA model parameters and an external estimate is unnecessary. Thus, for the case where a random walk model for the underlying process is valid, the appropriate level of smoothing is objectively determined.

As a demonstration, Pennington (1985) applied his noise reduction algorithm to groundfish trawl survey data for haddock (Melanogrammus aeglefinus) from the northeastern Atlantic coast of the United States. He found that the variances of the smoothed indices were "considerably lower" than those of the originals. However, this demonstration used an ARIMA model derived from a much longer time series that had been generated from a stock assessment based on commercial catch data. Pennington (1985) assumed that this model represented the underlying population and therefore did not develop models based on the observed time series. Although this assumption was perfectly reasonable, given that such alternative data (the stock assessment) were available, it cannot be applied to situations when only survey data are available to fishery analysts.

The ARIMA model Pennington (1985) derived from stock assessment results was a random walk model; therefore the appropriate level of noise reduction for the corresponding survey data could be objectively determined from the model parameters. Pennington's (1985) method was later used to apply random walk models to survey data (Fogarty et al. (1); Pennington, 1986; Anonymous, 1988, 1993). Pennington (1986) found that random walk models were appropriate for the survey time series considered in his study. However, random walk models were assumed a priori in the remaining three references (Fogarty et al. (1); Anonymous, 1988, 1993) to generate smoothed abundance trajectories; because less than 25 observations for each time series were considered in these references, reliable identification of the model structure for each time series was considered problematic and random walk models were used as "null" models.

When it is an appropriate description of the underlying process, a random walk model yields an objective determination of the degree of noise reduction appropriate to an observed time series. However, an a priori adoption of this model should be viewed with some skepticism. Additionally, if a random walk model is not an appropriate description of the underlying process, the resulting smoothed time series may seem reasonable, but the result no longer has support as the unobserved, underlying process. In this circumstance, we regard the effect of the ARIMA algorithm as merely smoothing, and not necessarily as noise reducing.

As such, we feel that the utility of ARIMA-based approaches to noise reduction for abundance indices derived from survey data has not been adequately explored to date. In addition, substantially longer time series (e.g., 40 observations) are now available with which to test this concept. In our study, we test the utility of the ARIMA time series noise reduction approach propounded by Pennington (1985), using time series of abundance indices from fishery-independent trawl survey data for nine finfish species (Table 1) during two seasons on Georges Bank. We first review the original methods developed by Cleveland and Tiao (1976) and Box et al. (1978). Framing the problem in terms of power spectra, we also offer some additional new insights into this noise reduction approach. Next, we apply the ARIMA-based noise reduction approach to the time series data and present the results. We have implemented Box et al.'s (1978) algorithm, not Pennington's (1985). Finally, we discuss our perceptions of the utility of this approach in light of our results and overall experience with it.

Materials and methods

General characteristics of ARIMA models

In this section, we first briefly review ARIMA models for stochastic processes. Then we review the approach of Box et al. (1978) for obtaining maximum likelihood estimates for an underlying ARIMA time series from a time series of observations with independent and identically distributed (IID) observation noise.

ARIMA models are parsimonious models that can adequately represent many stochastic time series (Box and Jenkins, 1976). Stochastic time series that can be represented by ARIMA models are essentially the output of a linear filter applied to an input time series of white noise (Box and Jenkins, 1976). We will refer to such time series as ARIMA processes.

For a zero-mean stochastic time series {[z.sub.t]} that can be expressed as an ARIMA model, we denote the model (using the notation of Box and Jenkins, 1976) as

[phi](B)[z.sub.t] = [alpha](B)[a.sub.t], (1)

where [z.sub.t] = the value of the time series at time t;

[phi](B) = the generalized autoregressive (AR) operator;

[alpha](B) = the moving average (MA) operator;

B = the backward shift operator; and

[a.sub.t] =IID normally distributed random variables with mean zero and variance [[sigma].sub.a.sup.2].

The backward shift operator B has the property that B [z.sub.t] = [z.sub.t-1]; hence [B.sup.m] [z.sub.t] = [z.sub.t-m]. The operators [phi](B) and [alpha](B) are polynomials in B of order p+d and q (respectively) such that

[phi](B) = 1 - [p+d.summation over (j=1)] [[phi].sub.j][B.sup.j]

[alpha](B) = 1 - [q.summation over (j=1)] [[alpha].sub.j][B.sup.j]. (2)

Furthermore, [phi](B) can be factored such that [phi](B)= [phi](B)[(1-B).sup.d], where the (ordinary) AR operator [phi](B) has a form similar to [phi](B) but is of order p. The operator [(1-B).sup.d] represents d sequential applications of the backward difference operator 1-B such that the original time series {[z.sub.t]} is self-differenced d times before application of the AR operator [phi](B).

As a shorthand, an ARIMA model that consists of an AR operator of order p, d backward difference operations, and a MA operator of order q is abbreviated as (p,d,q). A (p,0,0) model is referred to as an AR model of order p, or AR(p) in shorthand notation, and a (0,0,q) model is referred to as a MA model of order q (MA(q) in shorthand). A (p,0,q) model is referred to as an autoregressive-moving average model (ARMA(p,q) in shorthand) and a (0,d,q) model is referred to as an integrated moving average model (IMA(d,q) in shorthand). Finally, a random walk model is (0,1,0), while a random-walk-plus-uncorrelated-observation noise (RWPUN) model is (0,1,1).

In general, ARIMA models represent nonstationary time series. However, the time series that results from applying the backward difference operator d times to a (p,d,q) ARIMA process is a stationary ARMA(p,q) process. Typical constraints imposed on ARIMA models are that, when B is regarded as a complex variable, the polynomial in B representing the generalized autoregressive operator has zeros on or outside the unit circle ([absolute value of B]=1), whereas the polynomials representing the ordinary autoregressive operator and the moving average operator have zeros outside the unit circle. Additionally, the white noise variance, [[sigma].sup.2.sub.a], must be finite. These constraints ensure that a model is invertible and stationary (or, if not stationary, that at least it is of a suitable form).

One final property of ARIMA models is of interest here. The power spectrum for an ARMA (p,q) process represented by equation 1 is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (3)

where B = [e.sup.i2[pi]f; 0 [less than or equal to] f [less than or equal to] 1/2,

and the forward shift operator, F = [B.sup-1], has the opposite effect to B (i.e., [F.sup.m] [z.sub.t] = [z.sub.t+m]).

ARIMA models of time series with observation noise

Suppose, then, that [y.sub.1], [y.sub.2], ..., [y.sub.T] represent a time series of observations at times t = 1,2,... T such that

[y.sub.t] = [z.sub.t] + [e.sub.t], (4)

where the [z.sub.t]'s represent the unobserved, underlying process and the [e.sub.t]'s are IID normal variables with variance [[sigma].sub.e.sup.2] (i.e., e~N(0, [[sigma].sub.e.sup.2])) that are also independent of the [z.sub.t]'s (i.e., < [e.sub.j] [z.sub.k]>=0 for all j,k). The goal is to estimate the unobserved time series {[z.sub.t]} by using the observed time series {[y.sub.t]}. For the analysis of fishery-independent time series, it seems reasonable to assume that only the ARIMA model for the observed time series {[y.sub.t]} is known (it can be estimated using standard techniques). In particular, this assumption means that the model for {[z.sub.t]} is unknown within constraints implied by the observation equation. However, to develop the approach used here it is helpful to start as though the model for the unobserved process {[z.sub.t]} were known.

Thus, we assume that the time series {[z.sub.t]} can be represented by an ARIMA (p,d,q) process:

[phi](B)[z.sub.t] = [alpha](B)[c.sub.t], (5)

where the [c.sub.t]'s are IID and [c.sub.t] ~ N(0, [[sigma].sup.2.sub.c]). Substituting in [y.sub.t]-[e.sub.t] for [z.sub.t] and rearranging, one obtains

[phi](B)[y.sub.t] = [alpha](B)[c.sub.t] + [phi](B)[e.sub.t], (6)

which can be expressed as an ARIMA model for {[y.sub.t]} of the form

[phi](B)[y.sub.t] = [eta](B)[d.sub.t], (7)

where the [d.sub.t]'s are IID, [d.sub.t]~N(0, [[sigma].sup.2.sub.d]) and the MA operator is [eta](B). Thus, the generalized AR operator for the observed {[y.sub.t]} is identical to that for the unobserved {[z.sub.t]}. Furthermore, because [alpha](B) has order q and [phi](B) has order p+d, [eta](B) must be of order max(p+d,q). This requirement for [eta](B) constrains the order of potential ARIMA models that could describe the observed process: if p+d[less than or equal to]q then the observed process is also a (p,d,q) process, otherwise it is a (p,d,p+d) process.

In the more realistic situation where {[y.sub.t]} is observed, we can determine the generalized AR operator [phi](B) and MA operator [eta](B) for the observed process. Its ARIMA model will be order (P,D,Q), say, where the minimum possible value for Q is P+D. The model for the corresponding unobserved, underlying process {[z.sub.t]} will have order (p=P, d=D, q[less than or equal to]max(P+D,Q)) and its associated generalized AR operator will also be [phi](B). Furthermore, recognizing that [c.sub.t] and [e.sub.t] are independently distributed, it can be shown that the moving average operator for the unobserved process, [alpha](B), is additionally constrained by the ARIMA model for the observed process such that the following relationship must be satisfied:

[[sigma].sup.2.sub.c][alpha](B)[alpha](F) = [[sigma].sup.2.sub.d][eta](B)[eta](F) - [[sigma].sup.2.sub.e][phi](B)[phi](F). (8)

In this equation, [[sigma].sub.e.sup.2], [eta](B),and [phi](B) are known from the ARIMA model for the observed process, whereas [[sigma].sub.e.sup.2], [[sigma].sub.e.sup.2], and [alpha](B) are unknown.

In general, many combinations of [[sigma].sub.e.sup.2], [[sigma].sub.e.sup.2], and [alpha](B) will satisfy the equality. Defining an "acceptable model" for the unobserved process as one that, given the model for the observed process, [alpha]B) satisfies the previous equation and its zeros are on or outside the unit circle, Box et al. (1978) show that 1) for every given model of an observed process, at least one acceptable model for the unobserved process exists; 2) for a given model of an observed process, the possible values of [[sigma].sub.e.sup.2] are bounded; and 3) for a given model, every [[sigma].sub.e.sup.2] between 0 and the upper bound ([K.sup.*], say) determines a unique acceptable model. The upper bound on the observation error variance, [K.sup.*], is determined from the constraint that, for a model of the unobserved process to be acceptable, [[sigma].sub.c.sup.2][alpha](B)[alpha](F)[greater than or equal to] 0 everywhere on the unit circle (i.e., the power spectrum of the corresponding MA process is non-negative definite). Then, from equation 8, [K.sup.*] is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (9)

and is completely determined by the ARIMA model for the observed process. When [[sigma].sub.e.sup.2] = [K.sup.*], the variance of the added white noise is maximal, as will be the smoothing of the observed time series.

It is instructive to interpret Equations 8 and 9 in terms of constraints on the power spectra of {[z.sub.t]}, {[y.sub.t]}, and {[e.sub.t]}, although this interpretation is strictly correct only when {[z.sub.t]} and {[y.sub.t]} are stationary. Let [p.sub.z](f), [p.sub.y](f), [p.sub.e](f) denote the power spectra for {[z.sub.t]}, {[y.sub.t]}, and {[e.sub.t]}, respectively. Recalling the definition of the power spectrum (Eq. 3), Equation 8 on the unit circle can be easily recast (multiply both sides by 2/{[phi](B) [phi](F)}) as

[p.sub.z](f) = [p.sub.y](f) - [p.sub.e](f) (10)

because the power spectrum for white noise is constant with frequency: [p.sub.e](F)=2[[sigma].sup.2.sub.e]. Because power spectra are nonnegative definite and [p.sub.e](f) does not depend on f, the maximum possible observation noise variance [K.sup.*] corresponds to the minimum of [p.sub.y](f) over f (see Fig. 1). Thus, Equation 9 can be recast as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (11)

[FIGURE 1 OMITTED]

Note also that Equation 10 can be recast again as

[p.sub.z](f) = (1 - [p.sub.e](f) / [p.sub.y](f))[p.sub.y](f) = [[omega].sup.*](f)[p.sub.y](f), (12)

where [[omega].sup.*](f) represents the term in parentheses. It will be seen below that [[omega].sup.*](f) is identical to the polynomial of smoothing weights (Eq. 14) on the unit circle.

Returning to the estimation problem now, knowledge of [[sigma].sup.2.sub.e], then, together with the ARIMA model for the observed process, is sufficient to estimate {[z.sub.t]} (Box et al., 1978). When t is not close to the endpoints of the time series, the smoothed (maximum likelihood) estimate of [z.sub.t], which we denote as E([z.sub.t]|y), is a symmetric moving average filter of {[y.sub.t]} (Cleveland and Tiao, 1976; Box et al., 1978):

E([z.sub.t]|y) = [omega](B)[y.sub.t], (13)

such that

[omega](B = [[sigma].sup.2.sub.a]/[[sigma].sup.2.sub.d] [alpha](B)[alpha](F)/[eta](B)[eta](F) = [1 - [[sigma].sup.2.sub.e]/[[sigma].sup.2.sub.d] [phi](B)[phi](F)/[eta](B)[eta](F)], (14)

where the second equality follows from equation 8. Note that [omega](B) is identical to [[omega].sup.*](f) (Eq. 12) on the unit circle (B=[e.sup.-i2[pi]f)]), so that Equations 13 and 14 are also interpretable in terms of relations among power spectra. Equation 13 is equivalent to Equation 2 in Pennington (1985). It turns out (Cleveland and Tiao, 1976) that equation 14 is also valid for t near the endpoints of the time series; thus, one merely uses the ARIMA model for the observed time series to hindcast or forecast additional "observations" as needed. Box et al. (1978) describe a method for calculating the coefficients of w(B); because this reference may be difficult to obtain, we repeat their description in the Appendix.

In the case where the model for the underlying process is (0,1,0) (i.e., a random walk model), the model for the observed process is (0,1,1) (i.e., a RWPUN, model). For this case, Pennington (1985) noted that the value of the MA parameter [[eta].sub.1] of the observed process ([theta] in his notation) is

[[eta].sub.1] = [[sigma].sup.2.sub.e] / [[sigma].sup.2.sub.d], (15)

so that Equation 14 is completely determined by the ARIMA model for the observed process and the observation error variance, [[sigma].sup.2.sub.e], can be estimated.

Application of the noise reduction algorithm to bottom trawl survey data

We applied Box et al.'s (1978) noise reduction algorithm to 18 time series of abundance indices for finfish (nine speciesxtwo seasons; species are listed in Table 1) on Georges Bank in the northwest Atlantic. Time series for the fall survey spanned 40 years (1963-2002), and the spring time series spanned 36 years (1968-2003).

Stratified random bottom trawl surveys have been conducted annually on the northeastern continental shelf of the United States from Cape Hatteras to the Gulf of Maine in the fall since 1963 and in the spring since 1968 by the National Marine Fisheries Service (NMFS), Northeast Fisheries Science Center (NEFSC) (Azarovitz, 1981; Anonymous, 1988; Reid et al., 1999). We derived annual time series of estimated total abundance during the spring and fall seasons for nine finfish species from trawl survey catch records for the Georges Bank area. Reported catches (biomass) from strata encompassing Georges Bank were expanded on a per-tow basis by using species and length-specific corrections for catchability and coefficients from Edwards (1968) and Harley and Myers (2001). Annual stratified mean (expanded) catch-per-tow was then calculated for each species for both seasonal surveys. Finally, annual total population abundance was calculated by applying a swept area factor to the stratified mean catch-per-tow. Abundance indices corresponding to fall surveys spanned 1963-2002, and those corresponding to spring surveys spanned 1968-2003. Following Pennington (1985), we assumed that changes in catchability induced a lognormal error structure on the observed time series. However, because some time series (notably those for the schooling pelagic fish herring and mackerel) included zeroes, a ln(x) transformation was not applicable to all series. Thus, before further analysis, all time series were ln(x+1) transformed.

For each resulting time series, we used SAS (vers. 8, SAS Institute, Cary, NC) to identify candidate ARIMA models, to estimate parameters for these models, and to extend the time series by using forecasts and hindcasts before applying the smoothing algorithm. Candidate model structures for each time series were based on examination of empirical autocorrelation, partial autocorrelation, and inverse autocorrelation functions for the series (Box and Jenkins, 1976). When these functions indicated that the series was nonstationary, we applied the backward difference operator (1-B) to the series and examined the correlation functions for the new (differenced) series. In addition, because of its special significance in terms of interpretation, we always tried a RWUPN model as a candidate. Once candidate models were identified, we estimated coefficients for each model and calculated the associated Akaike information criterion (AIC; Akaike, 1973). AIC provides an objective criterion based on information theory for selecting the "best" approximating model from among a group of good candidate models (i.e., the criterion selects the model with the minimum AIC).

We examined the residuals and the empirical autocorrelation, inverse autocorrelation, and partial autocorrelation functions of the residual time series for significant deviation from white noise. When significant deviation was indicated, we dropped the model from further consideration. We also dropped models with orders (P,D,Q) that were inconsistent with the assumption of additive (after log transformation) observation noise (Q[greater than or equal to]P+D). Following this screening procedure, we were left with a group of "good" alternative models. We then selected the ARIMA model with the smallest AIC from among the remaining candidates as the "best" model to smooth the data.

We developed a MATLAB (vers. 6.5, The Mathworks Inc., Natick, MA) program to perform the noise reduction for each time series based on an extended version of the series and its associated ARIMA model. To allow smoothed estimates to be calculated near the ends of each time series, we extended each time series before smoothing to 40 years before its start by hindcasting with the selected ARIMA model and to 40 years past its end by forecasting with the model. Using the MATLAB program, we then calculated [K.sup.*], the maximal smoothing weights [omega](B) corresponding to [[sigma].sub.e.sup.2] = [K.sup.*], and the smoothed estimates E([z.sub.t]|y), following the approach of Box et al. (1978) outlined previously.

Results

The abundance indices we derived from fishery-independent bottom trawl surveys for nine finfish species during two seasons displayed a wide variety of trends, as well as a high degree of apparent interannual variability. This variability may be associated with environmentally driven, high-frequency changes in catchability, but we regarded it here as observation noise (Figs. 2 and 3). For example, springtime winter skate (Leucoraja ocellata) biomass appeared to increase by a factor of six during the early 1980s from a lower (but highly variable) mean level of ~140,000 t, to which it subsequently returned in the early 1990s (Fig. 2A). Yellowtail flounder (Limanda ferruginea), in contrast, declined by a factor of 10 during the 1970s and early 1980s from a high at the beginning of the time series and then began to rebound in the latter 1980s (Fig. 2G). Most recently, yellowtail flounder appear close to regaining their earlier peak abundance.

[FIGURES 2-3 OMITTED]

The 18 ARIMA models corresponding to the ln(x+1)-transformed abundance indices formed a diverse set (Tables 2 and 3). Half the time series were found to be nonstationary; however, one application of backward differencing (d=1) sufficed to achieve stationarity in all nine instances. Interestingly, the ARIMA models for all nine of these time series were consistent with RWPUN models (i.e., [0,1,1] models). For all the models, autoregressive orders (p) ranged between 0 and 3, and moving average orders (q) ranged from 1 to 5. The spring and fall models for little skate (L. erinacea) had the highest AR order, and the spring model for silver hake (Merluccius bilinearis) had the highest MA order. Only the models for little skate, Atlantic herring (Clupea harengus), and yellowtail flounder exhibited the same ARIMA order for both the spring and fall. And although none of the models reflected pure AR processes (q=0, an impossibility given the observation noise assumption), two were found to be pure MA processes (winter flounder and Atlantic mackerel [Scomber scombrus], both in spring). All nine IMA processes were RWPUN processes (i.e., (0,1,1)).

The effect of ARIMA-based noise reduction on these 18 sets of indices was fairly variable in outcome (Figs. 2 and 3). Little or no smoothing occurred for little skate (spring, Fig. 2B), silver hake (spring, Fig. 2D), and haddock (Melanogrammus aeglefinus) (fall, Fig. 3F). Moderate smoothing occurred for winter skate (spring, Fig. 2A; fall, Fig. 3A), little skate (fall, Fig. 3B), Atlantic herring (fall, Fig. 3C), haddock (spring, Fig. 2F), yellowtail flounder (spring, Fig. 2G; fall, Fig. 3G), winter flounder (Pseudopleuronectes americanus) (fall, Fig. 3H), and Atlantic mackerel (spring, Fig. 2I). A high degree of smoothing occurred for Atlantic herring (spring, Fig. 2C), silver hake (fall, Fig. 3D), Atlantic cod (Gadus morhua) (spring, Fig. 2E; fall, Fig. 3E), winter flounder (spring, Fig. 2H), and Atlantic mackerel (fall, Fig. 3I).

Impressions obtained from visually comparing the original and smoothed time series were consistent with the fraction of innovation variance ([K.sup.*]/[[sigma].sub.d.sup.2], Tables 2 and 3) in the observed time series ascribed to observation error. Recall that we chose to perform the maximal amount of smoothing possible (i.e., taking [[sigma].sub.d.sup.2]=[K.sup.*] in Eq. 13). Our visual classifications fell fairly neatly into the following categories based on [K.sup.*]/[[sigma].sub.d.sup.2]: 1) a low degree of smoothing corresponded to [K.sup.*]/[[sigma].sub.d.sup.2]<0.25, 2) a moderate degree of smoothing corresponded to 0.25< [K.sup.*]/[[sigma].sub.d.sup.2]<0.60, and 3) a high degree of smoothing corresponded to 0.60[greater than or equal to][K.sup.*]/[[sigma].sub.d.sup.2].

In addition, it appeared that the degree of smoothing varied inversely with the order of the moving average component of the ARIMA model for the observed time series. When the order of the moving average component was greater than 2, the degree of overall smoothing was typically small. Substantially more smoothing was evident when the order of the moving average component was 1, because of the equivalence with an exponentially weighted smoother.

The nine time series found to be RWPUN processes provided a means to check whether our choice to perform the maximal amount of smoothing was reasonable. For such time series, the value of the moving average coefficient, [[eta].sub.1], is equal to the ratio of the observation noise variance, [[sigma].sub.e.sup.2], to the innovation variance [[sigma].sub.d.sup.2]. Consequently, [[sigma].sup.2.sub.d]x[[eta].sub.1]=[[sigma].sup.2.sub.e]. For eight of the nine series, the ratio of [[sigma].sub.d.sup.2] to [K.sup.*] was greater than 80% (Table 4).

Discussion

Abundance indices derived from large-scale fishery-independent surveys typically exhibit interannual variability much higher than one would expect from within-survey variance (Pennington, 1985). Although true variability in the underlying population due to population-dynamic processes is reflected in the variability of an index, so too is observational noise arising both from within-survey sampling variability as well as from environmentally driven factors that affect catchability. Low signal-to-noise ratios in abundance indices due to high observational noise reduce one's ability to detect important changes or trends in actual population abundance.

To reduce the impact of white observation noise on time series data, Cleveland and Tiao (1976) developed an approach to noise reduction and smoothing that was based on their knowledge of the ARIMA model (Box and Jenkins, 1976) for an associated unobserved but underlying stochastic process. Box et al. (1978) extended this approach to address the situation where the ARIMA model for the underlying process was unknown, relying instead on an ARIMA model associated with the observed time series. In general, and certainly in regard to fishery-independent survey data, a model structure for the unobserved, underlying process will not be available. Hence, Box et al.'s (1978) approach will be the norm.

In the situation where the observed time series is stationary, we found that a frequency domain interpretation of Box et al.'s (1978) algorithm is particularly enlightening. When the times series is stationary, observation noise increases the power spectral density (PSD) of the observed process over that of the unobserved process by a fixed amount at all frequencies (Fig. 1). Consequently, the PSD of the unobserved process has the same shape as the PSD of the observed process, but with a fixed amount removed at all frequencies. Because the unobserved PSD must be nonnegative, the maximum observation noise (and thus the maximum noise reduction and smoothing that may occur) consistent with additive white noise corresponds to the minimum value of the observed process's PSD. Presumed higher values for the observation noise would require that the PSD for the unobserved process be negative over some frequency range. Box et al.'s (1978) algorithm computes the maximum possible observation noise and uses a time domain formula to calculate a smoothed, "unobserved" time series consistent with additive white noise for any noise level up to the maximum.

The ARIMA-based noise reduction approach was first applied to fisheries trawl survey data by Pennington (1985), who developed an alternative algorithm to that of Box et al. (1978) for estimating the smoothed time series. This algorithm was subsequently used in several studies (Pennington, 1985; Fogarty et al. (1); Pennington, 1986; Anonymous, 1988, 1993) to smooth time series of abundance indices derived from trawl survey data for the northeast coast of the United States. Of these studies, perhaps only Pennington (1986) constitutes a convincing demonstration of the utility of the ARIMA-based approach to time series noise reduction when the ARIMA model for the underlying process is unknown. Pennington's (1985) demonstration relied on an ARIMA model developed for the (usually unobserved) underlying population that was based on stock assessment results for the particular species considered. In contrast, RWPUN models were assumed a priori in Fogarty et al. (1) and Anonymous (1988, 1993) because the short time series (<25 observations per series) considered made identification of the underlying ARIMA models problematic. Only in Pennington (1986) was a model used that was fitted to trawl survey data and tested for appropriateness.

However, substantially longer time series are now available to test the ARIMA-based noise reduction concept. Consequently, to revisit the utility of the ARIMA-based approach for smoothing time series data derived from fisheries-independent trawl survey data, we applied Box et al.'s (1978) approach to smoothing to time series of annual abundance indices derived from the NMFS/NEFSC fisheries-independent bottom trawl survey of nine finfish species during two seasons on Georges Bank. Time series for the fall spanned 40 years (1963-2002), and the spring time series spanned 36 years (1968-2003). The species, comprising two elasmobranchs, three groundfish, two flatfish, and two pelagic schooling species, presented a broad range of life history characteristics (Table 1).

The noise reduction results we obtained varied among species and between seasons. Despite smoothing at the maximum level of noise reduction consistent with each model, very little smoothing was obtained for haddock (fall), little skate (spring), and silver hake (spring). The models for these time series had among the highest moving average orders (3-5). Examination of the PSD for each model in the frequency domain revealed very small minima, indicating no scope for noise reduction. Typically, models that had a MA order >2 exhibited substantially less smoothing than models with a MA order [less than or equal to]2. Models that were of MA order 1 generally resulted in the greatest smoothing. Models that were of MA order 0 did not (and could not) occur.

Of the 18 time series we considered (Tables 2 and 3, Figs. 2 and 3), only half were adequately represented as random-walk-plus-uncorrelated-noise (RWPUN) models. The ARIMA models we developed were varied in structure, ranging from a simple MA(1) model to rather complicated models with multiple parameters. Thus, our results provide evidence against the appropriateness of assuming a particular model structure a priori when the objective of the analysis is to identify the underlying dynamic structure of the population. This evidence is further strengthened by the results of Becerra-Munoz et al. (1999), who found only 9 of 52 abundance time series for finfish species from the NMFS/NEFSC bottom trawl survey that corresponded to random walk models.

As an exercise, we also attempted to smooth the nine data sets that were not adequately described as RWPUN models, using this model structure as an a priori assumption, even though our analysis indicated that other models were more appropriate. We were not able to estimate convergent models for three species: Atlantic cod (fall), winter flounder (spring), and Atlantic mackerel (spring). For the remaining six time series, the smoothed results appeared to be quite reasonable (Fig. 4), although we obtained little noise reduction when we employed the "correct" ARIMA model. The RWPUN-smoothed time series for haddock (Fig. 4, C and F) were similar to that for spawning biomass derived from virtual population analysis (see Brodziak et al. (2)), but the smoothed time series for silver hake (Fig. 4E) exhibited higher frequency variability than that found for total biomass with a production model (see Brodziak et al. (3)). From the standpoint of estimating the unobserved underlying process, these smoothed results should be viewed with some skepticism: the use of the RWPUN model is rather arbitrary in this situation and it may impose artificial structure on the smoothed results. However, it may be that these time series do not meet one of the key assumptions of the noise reduction method: namely that the observation noise is uncorrelated. The ARIMA models for all six time series had MA orders [greater than or equal to]3, and one effect of correlated observation noise could be to increase the MA order of the observed time series beyond that expected for uncorrelated noise.

[FIGURE 4 OMITTED]

Of the nine species considered, only three (little skate, Atlantic herring, and flounder) had models that exhibited the same ARIMA order for both the fall and spring surveys. Taken at face value, this would indicate that the other six species exhibited substantially different dynamical processes during the fall and spring that influenced abundance on Georges Bank. One potential mechanism for this could be differential seasonal migration patterns that result in changes in catchability that have different autocorrelation structures. For example, cod are distributed across the bank during the spring survey, but are found only in deeper waters on the periphery of the bank during the fall where they may be less available to the survey. Assuming that all Georges Bank cod are available to the spring survey, if the fraction of cod available to the survey in the fall is density dependent or is driven by autocorrelated environmental conditions, then the fall survey abundance will exhibit dynamical behavior different from that of the spring survey abundance (note: if this were the case, it would be inappropriate to smooth the fall time series for cod with the ARIMA noise reduction approach applied in our study because, as noted previously, one of the basic assumptions with this approach is that the observation noise is uncorrelated). Other plausible mechanisms can be developed, as well. However, we feel it more likely that the inconsistency in ARIMA order between spring and fall surveys for the same species is an estimation problem and indicates that even a 40-year time series may not be long enough to reduce the variability inherent in ARIMA model estimation to reasonable levels for most species.

Alternative methods, such as locally weighted scatterplot smoothing (LOESS), moving average filters, exponential smoothing filters, Kalman filters, and frequency-domain approaches can be applied to time series to achieve smoother results (e.g., Cleveland and Grosse, 1991; Hamilton, 1994). These approaches typically employ at least one user-determined parameter that can be used to change the amount of smoothing that an algorithm achieves. Generally, one "fiddles" with the adjustable parameters until a "nice," smoothed fit is achieved. However, we think it important to distinguish between these smoothing algorithms and the ARIMA-based noise reduction algorithms. It is quite possible to smooth out real fluctuations in the underlying population process. The principal advantage that we see for the ARIMA-based noise reduction algorithms (used with an appropriate model) over alternative methods is that the former provide a more objective approach to determining an appropriate level of smoothing. As noted previously, Pennington (1985) showed that when a RWPUN model is appropriate, the ARIMA smoothing approach is completely determined by the ARIMA model for the observed time series because it is possible to determine the observation noise variance from the model parameters. In the more general case, Box et al.'s (1978) algorithm at least yields a maximum value for the variance of the observational noise and thus sets an upper limit to the amount of noise reduction and smoothing that can be achieved. For trawl survey data, our results from nine time series where RWPUN models were appropriate (and we can consequently estimate the actual observation noise variance) indicate that smoothing at ~90% of the maximum possible noise reduction level is not an unreasonable default percentage (Table 4).

One drawback to the greater application of ARIMA-based noise reduction methods to time series data is the lack of an integrated software package that allows a user 1) to quickly evaluate an appropriate ARIMA model for a given time series, and 2) to calculate the smoothed time series. We used SAS for the first step and MATLAB for the second, but we found this arrangement rather awkward and burdensome. However, econometrically oriented software packages such as ForecastPro (4) or AutoBox (5) that automate model selection may substantially simplify the first step even if they don't address the second step.

On the whole, ARIMA-based time series models appear to provide the basis for a more objective approach to reducing observation noise in time series data, including time series of fishery abundance indices derived from trawl survey data, than do more conventional smoothing approaches. In the absence of additional information regarding the level of observation noise, we recommend smoothing trawl survey data at 90% of the maximum possible noise reduction level. We also suggest that development of an integrated software package for implementing ARIMA-based noise reduction will facilitate future use of this method.

Finally, if a smoothed time series is desired (e.g., for graphical presentation only), then use of a RWPUN model in lieu of a model-fitting exercise will generally yield a curve pleasing to the eye. Alternatively, other methods such as LOESS could be employed to generate the smoothed results. However, if the resulting time series is to be used for further analysis of the dynamical behavior of the fish stock, we strongly recommend that a model-fitting approach be used to identify the most appropriate ARIMA model for the observed time series, from which the time series for the unobserved, underlying process can be computed. Otherwise, real fluctuations in the underlying process may be oversmoothed, resulting in an apparent dynamical behavior that displays little variability. This oversmoothing, in turn, may lead to erroneous conclusions being drawn regarding, for example, the resiliency of a stock to exploitation or environmental change, and to perhaps concomitant errors being propagated in advice provided to fishery managers.

Appendix

For the convenience of the reader, we summarize here Box et al.'s (1978) algorithm to calculate the coefficients of the smoothing polynomial [omega](B). Recall from Equation 14 that

[omega](B) = [1 - [[sigma].sup.2.sub.e]/[[sigma].sup.2.sub.d] [phi](B)[phi](F)/[eta](B)[eta](F], (A.1)

where [phi](B) is a polynomial of order P+D and [eta](B) is a polynomial of order Q. Because Q [greater than or equal to] P+D, one can write

[phi](B) = 1 - [[phi].sub.1][B.sup.1] - [[phi].sub.2][B.sup.2] ... - [[phi].sub.Q][B.sup.Q]; [eta](B) = 1 - [[eta].sub.1][B.sup.1] - [[eta].sub.2][B.sup.2] ... - [[eta].sub.Q][B.sup.Q], (A.2)

where [[phi].sub.j] = 0 for j>P+D.

First, define

C(B) [equivalent to] [phi](B)/[eta](B); C(B) = 1 + [C.sub.1][B.sup.1] + [C.sub.2][B.sup.2] +.... (A.3)

One can solve for the coefficients of C using an iterative process by recognizing that the coefficients of each power of B in the following expression must be zero:

0 = [phi](B)- [eta](B)C(B).

Consequently, one obtains

[C.sub.1] = [[eta].sub.1] - [[phi].sub.1]; [C.sub.j] = [[eta].sub.j] - [[phi].sub.j] + [j-1.summation over (i=1)] [C.sub.j-i][[eta].sub.i], j = 2,... Q (A.4)

and higher order coefficients of C (i.e., for j>Q) can be computed recursively using the relation

[C.sub.j] = [Q.summation over (i=1)] [C.sub.j-i][[eta].sub.i]. (A.5)

Next, define

X(B,F) [equivalent to]C(B)[[phi](F)/[eta](F) (A.6)

where X(B,F) = [X.sub.0] + [[infinity].summation over (j=1)] [X.sub.j] ([B.sup.j] + [F.sup.j]).

From this definition, one obtains

X(B, F)[eta](F) [equivalent to]C(B)[phi](F). (A.7)

The largest degree of F on the righthand side of this equation is Q. By matching coefficients of [F.sup.j] on both sides of the equation, Box et al. (1978) found, for j=0,1,... Q, that

[eta]X = C[phi], (A.8)

where X, [phi] are Q+1 column vectors [X.sup.T]=([X.sub.0], [X.sub.1],..., [X.sub.Q]), [[phi].sup.T]=(1, -[[phi].sub.1],..., -[[phi].sub.Q]) and [eta], C are (Q+1)x(Q+1) matrices such that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A.9)

Thus, Equations A.8 and A.9 define a set of Q+1 equations in Q+1 unknowns ([x.sub.0], [X.sub.1],..., [X.sub.Q]) which can be solved for by matrix inversion such that

X = [n.sup.-1]C[phi]. (A.10)

For j>Q, the [X.sub.j]'s are computed recursively by using the relation

[X.sub.j] = [Q.summation over (i=1)] [X.sub.j-1][[eta].sub.i]. (A.11)

Finally, the coefficients of [omega](B) are given by

[[omega].sub.0] = 1 - [[sigma].sup.2.sub.e]/[[sigma].sup.2.sub.d] [X.sub.0]; [[omega].sub.j] = -[[sigma].sup.2.sub.e]/[[sigma].sup.2.sub.d] [X.sub.j], j = 1,.... (A.12)

Acknowledgments

We would like to thank M. Pennington for his insightful comments on an early version of this manuscript. The comments and suggestions from two anonymous reviewers were also very helpful and are much appreciated. This work was supported by the National Research Council, which provided funding support for one of us (Stockhausen) as a postdoctoral fellow at the Northeast Fisheries Science Center.

Manuscript submitted 12 November 2004 to the Scientific Editor's Office.

Manuscript approved for publication 18 April 2006 by the Scientific Editor.

Literature cited

Akaike, H. 1973. Information theory as an extension of the maximum likelihood principal. In Second international symposium on information theory (B. N. Petrov, and F. Csaki, eds.), 267-281. Akademiai Kiado, Budapest, Hungary.

Anonymous. 1988. An evaluation of the bottom trawl survey program of the Northeast Fisheries Center. NOAA Tech. Memo. NMFS-F/NEC-52, 83 p. Northeast Fish. Center, 166 Water St., Woods Hole, MA 02543.

1993. Status of fishery resources off the northeastern United States for 1993. NOAA Tech. Memo. NMFSF/NEC-101, 140 p. Northeast Fish. Center, 166 Water St., Woods Hole, MA 02543.

Azarovitz, T. R. 1981. A brief historical review of the Woods Hole Laboratory trawl survey time series. In Bottom trawl surveys (W. G. Doubleday, and D. Rivard, eds.), 62-67. Can. Spec. Pub. Fish. Aquat. Sci., vol. 58.

Becerra-Munoz, S., D. B. Hayes, and W. W. Taylor. 1999. Stationarity and rate of damping of modeled indices of fish abundance in relation to their exploitation status in the Northwest Atlantic Ocean. Ecol. Mod. 117:225-238.

Box, G. E. P., S. C. Hillmer, and G. C. Tiao. 1978. Analysis of seasonal time series. In Seasonal analysis of economic time series (A. Zellner, ed.), 309-334. U.S. Dep. Comm., Washington, DC.

Box, G. E. P., and G. M. Jenkins. 1976. Time series analysis: forecasting and control, 575 p. Holden-Day, Oakland, CA.

Byrne, C. J., T. R. Azarovitz, and M. P. Sissenwine. 1981. Factors affecting variability of research trawl surveys. Can. Spec. Pub. Fish. Aquat. Sci. 58:238-273.

Clark, S. H. 1979. Application of bottom trawl survey data to fish stock assessment. Fisheries 4:9-15.

Cleveland, W. S., and E. Grosse. 1991. Computational methods for local regression. Stats. and Comp. 1:47-62.

Cleveland, W. P., and G. C. Tiao. 1976. Decomposition of season time series: a model for the Census X-11 program. J. Am. Stat. Assoc. 71:581-587.

Collie, J. S., and M. P. Sissenwine. 1983. Estimating population size from relative abundance data measured with error. Can. J. Fish. Aquat. Sci. 40:1971-1983.

Edwards, R. L. 1968. Fishery resources of the North Atlantic area. In The future of the fishing industry in the United States (D. W. Gilbert, ed.), p. 52-60. Univ. Washington, Seattle, WA.

Enders, W. 2004. Applied econometric time series, 460 p. John Wiley and Sons, Inc., Hoboken, NJ.

Grosslein, M. D. 1969. Groundfish survey program of BCF Woods Hole. Comm. Fish. Rev. 31:22-30.

Hamilton, J. 1994. Time series analysis, 820 p. Princeton Univ. Press, Princeton, NJ.

Harley, S. H., and R. A. Myers. 2001. Hierarchical Bayesian models of length-specific catchability of research trawl surveys. Can. J. Fish. Aquat. Sci. 58:1569-1584.

Helser, T. E., and D. B. Hayes. 1995. Providing quantitative management advice from stock abundance indices based on research surveys. Fish. Bull. 93:290-298.

Lloret, J., J. Lleonart, and I. Sole. 2000. Time series modeling of landings in northwest Mediterranean Sea. ICES J. Mar. Sci. 57:171-184.

Nogueira, E., F. F. Perez, and A. F. Rios. 1998. Modelling nutrients and chlorophyll a time series in an estuarine upwelling ecosystem (Ria de Vigo: NW Spain) using the Box-Jenkins approach. Estuar. Coast. Shelf Sci. 46:267-286.

Pennington, M. 1985. Estimating the relative abundance of fish from a series of trawl surveys. Biometrics 41:197-202.

1986. Some statistical techniques for estimating abundance indices from trawl surveys. Fish. Bull. 84: 519-525.

Reid, R. N., F. P. Almeida, and C. A. Zetlin. 1999. Fishery-independent surveys, data sources, and methods. NOAA Tech. Memo. NMFS-NE-122. Northeast Fish. Sci. Center, 166 Water St., Woods Hole, MA 02543.

(1) Fogarty, M. J., J. S. Idoine, F. P. Almeida, and M. Pennington. 1986. Modeling trends in abundance based on research vessel surveys. ICES CM (council meeting) 1986/G, p. 92. ICES, Copenhagen, Denmark.

(2) Brodziak, J., M. Traver and L. Col. 2005. Georges Bank haddock. In Assessment of 19 northeast groundfish stocks through 2004 (R. K. Mayo, and M. Terceiro, eds.), section 2, p. 30-80. 2005 groundfish assessment review meeting. Northeast Fisheries Science Center, Woods Hole, Massachusetts; 15-19 August 2005. NEFSC Ref. Doc. 05-13. NEFSC, 166 Water Street, Woods Hole, MA 02543.

(3) Brodziak, J. K. T., E. M. Holmes, K. A. Sosebee, and R. K. Mayo. 2001. Assessment of the silver hake resource in the Northwest Atlantic in 2000, 134 p. NEFSC Ref. Doc. 0103. NEFSC, 166 Water Street, Woods Hole, MA 02543.

(4) Business Forecast Systems, Inc. 2006. Website: http://www. forecastpro.com/products/fpfamily/index.html (accessed on 29 March 2006).

(5) Automatic Forecasting Systems. 2003. Website: http:// www.autobox.com/autoboxdesc.htm (accessed on 29 March 29 2006).

William T. Stockhausen (contact author)

Michael J. Fogarty

Email address for W. T. Stockhausen: William.Stockhausen@noaa.gov

National Ocean and Atmospheric Administration

National Marine Fisheries Service

Northeast Fisheries Science Center

166 Water Street

Woods Hole, MA 02543

Present address for corresponding author: National Marine Fisheries Service

Alaska Fisheries Science Center

7600 Sand Point Way NE

Seattle, Washington 98115

**********

Time series of species abundance from fishery-independent surveys, such as bottom trawl or acoustic surveys, are important in monitoring temporal change in the abundance of marine populations. For commercially important species, catch and effort data from the commercial fishery may be available, allowing estimation of temporal trends of the stock population by means of stock assessment models (e.g., virtual population analysis). However, such records are not available for many species, especially those with little commercial (but perhaps significant ecological) value. Fishery-independent surveys may thus constitute the only source of information for assessing temporal changes in the abundance of these species (Pennington, 1985; Helser and Hayes, 1995).

Annual estimates of abundance derived from fisheries-independent surveys are typically regarded as providing a relative measure of population abundance (i.e., they are indices of abundance, not true estimates of total population size) (Grosslein, 1969; Clark, 1979). Thus, the expected value of the abundance index (e.g., mean catch-per-tow for trawl surveys) is regarded as proportional to the size of the actual population, although the constant of proportionality (the catchability) is unknown. As such, relative changes in an abundance index should reflect similar relative changes in the actual population, and trends in the time series of such an index should reflect similar trends in the corresponding population.

Unfortunately, abundance indices derived from large-scale fishery-independent surveys typically exhibit interannual variability much higher than one would expect from withinsurvey variance (Byrne et al., 1981; Pennington, 1985). Part of the variability in such indices is presumably due to the variability in the underlying population--a variability that is caused by population-dynamic processes such as recruitment. However, part of the variability is due to observation noise that arises from both within-survey sampling variability because of the heterogeneous distribution of many fish stocks (Byrne et al., 1981), and because of environmentally driven factors that affect catchability over time (Byrne et al., 1981; Collie and Sissenwine, 1983). Low signal-to-noise ratios in abundance indices that are due to high observation noise reduce chances of detecting important changes or trends in actual population abundance. Variability due to within-survey sampling can be reduced (before the fact) by adding more stations to a survey, but additional stations will not reduce variability caused by changes in catchability.

Time series modeling using autoregressive integrated moving average (ARIMA; Box and Jenkins, 1976) models provides an approach to removing observation noise from abundance estimates. ARIMA models are frequently used in economic forecasting (Enders, 2004) and are becoming more common in fisheries research. Recent applications of ARIMA models to other fisheries problems include forecasting monthly landings in the Mediterranean (Lloret et al., 2000), testing theories of population dynamics (Becerra-Munoz et al., 1999), and modeling nutrient dynamics in an upwelling system (Nogueira et al., 1998).

In the context of reducing the influence of observational noise in time series data, Cleveland and Tiao (1976) first developed a noise-reduction and smoothing algorithm for processes that could be described by an ARIMA time series model. Their approach requires that the ARIMA model for the unobserved, underlying process be known. This known model, in turn, uniquely determines the ARIMA model for the observed time series contaminated by observation noise and allows one to estimate the variance of the observation noise. Unfortunately, although the ARIMA model for an unobserved, underlying process may be known in some instances (from theory, perhaps), in many cases the model for the unobserved process will be unknown.

Box et al. (1978) extended Cleveland and Tiao's (1976) ideas and developed a noise-reduction algorithm based on the ARIMA model for the observed time series. However, the ARIMA model for the observed time series merely constrains, but does not determine, the model for the unobserved, underlying process; it provides only an upper bound for the observation error variance. Consequently, this approach generally requires an external estimate of the observation error variance to determine the appropriate level of noise reduction.

Pennington (1985) first applied these ARIMA-based time series modeling techniques to smoothing abundance indices derived from trawl survey data. He assumed that an observed abundance time series reflected a combination of the underlying population abundance and independent, uncorrelated, and multiplicative observation noise (the latter arising perhaps from environmentally driven changes in catchability). He further assumed that both the (log-transformed) observed time series and unobserved population process could be represented by ARIMA models. Pennington (1985) then developed an alternative algorithm to that of Box et al. (1978); his derivation allowed particular simplification in the case where the underlying population process could be modeled as a random walk. In this simple case, the resulting noise reduction filter is an exponentially-weighted average of the observed time series for the endpoint of the time series (Pennington, 1985). More importantly, the observation error variance can be easily estimated from the ARIMA model parameters and an external estimate is unnecessary. Thus, for the case where a random walk model for the underlying process is valid, the appropriate level of smoothing is objectively determined.

As a demonstration, Pennington (1985) applied his noise reduction algorithm to groundfish trawl survey data for haddock (Melanogrammus aeglefinus) from the northeastern Atlantic coast of the United States. He found that the variances of the smoothed indices were "considerably lower" than those of the originals. However, this demonstration used an ARIMA model derived from a much longer time series that had been generated from a stock assessment based on commercial catch data. Pennington (1985) assumed that this model represented the underlying population and therefore did not develop models based on the observed time series. Although this assumption was perfectly reasonable, given that such alternative data (the stock assessment) were available, it cannot be applied to situations when only survey data are available to fishery analysts.

The ARIMA model Pennington (1985) derived from stock assessment results was a random walk model; therefore the appropriate level of noise reduction for the corresponding survey data could be objectively determined from the model parameters. Pennington's (1985) method was later used to apply random walk models to survey data (Fogarty et al. (1); Pennington, 1986; Anonymous, 1988, 1993). Pennington (1986) found that random walk models were appropriate for the survey time series considered in his study. However, random walk models were assumed a priori in the remaining three references (Fogarty et al. (1); Anonymous, 1988, 1993) to generate smoothed abundance trajectories; because less than 25 observations for each time series were considered in these references, reliable identification of the model structure for each time series was considered problematic and random walk models were used as "null" models.

When it is an appropriate description of the underlying process, a random walk model yields an objective determination of the degree of noise reduction appropriate to an observed time series. However, an a priori adoption of this model should be viewed with some skepticism. Additionally, if a random walk model is not an appropriate description of the underlying process, the resulting smoothed time series may seem reasonable, but the result no longer has support as the unobserved, underlying process. In this circumstance, we regard the effect of the ARIMA algorithm as merely smoothing, and not necessarily as noise reducing.

As such, we feel that the utility of ARIMA-based approaches to noise reduction for abundance indices derived from survey data has not been adequately explored to date. In addition, substantially longer time series (e.g., 40 observations) are now available with which to test this concept. In our study, we test the utility of the ARIMA time series noise reduction approach propounded by Pennington (1985), using time series of abundance indices from fishery-independent trawl survey data for nine finfish species (Table 1) during two seasons on Georges Bank. We first review the original methods developed by Cleveland and Tiao (1976) and Box et al. (1978). Framing the problem in terms of power spectra, we also offer some additional new insights into this noise reduction approach. Next, we apply the ARIMA-based noise reduction approach to the time series data and present the results. We have implemented Box et al.'s (1978) algorithm, not Pennington's (1985). Finally, we discuss our perceptions of the utility of this approach in light of our results and overall experience with it.

Materials and methods

General characteristics of ARIMA models

In this section, we first briefly review ARIMA models for stochastic processes. Then we review the approach of Box et al. (1978) for obtaining maximum likelihood estimates for an underlying ARIMA time series from a time series of observations with independent and identically distributed (IID) observation noise.

ARIMA models are parsimonious models that can adequately represent many stochastic time series (Box and Jenkins, 1976). Stochastic time series that can be represented by ARIMA models are essentially the output of a linear filter applied to an input time series of white noise (Box and Jenkins, 1976). We will refer to such time series as ARIMA processes.

For a zero-mean stochastic time series {[z.sub.t]} that can be expressed as an ARIMA model, we denote the model (using the notation of Box and Jenkins, 1976) as

[phi](B)[z.sub.t] = [alpha](B)[a.sub.t], (1)

where [z.sub.t] = the value of the time series at time t;

[phi](B) = the generalized autoregressive (AR) operator;

[alpha](B) = the moving average (MA) operator;

B = the backward shift operator; and

[a.sub.t] =IID normally distributed random variables with mean zero and variance [[sigma].sub.a.sup.2].

The backward shift operator B has the property that B [z.sub.t] = [z.sub.t-1]; hence [B.sup.m] [z.sub.t] = [z.sub.t-m]. The operators [phi](B) and [alpha](B) are polynomials in B of order p+d and q (respectively) such that

[phi](B) = 1 - [p+d.summation over (j=1)] [[phi].sub.j][B.sup.j]

[alpha](B) = 1 - [q.summation over (j=1)] [[alpha].sub.j][B.sup.j]. (2)

Furthermore, [phi](B) can be factored such that [phi](B)= [phi](B)[(1-B).sup.d], where the (ordinary) AR operator [phi](B) has a form similar to [phi](B) but is of order p. The operator [(1-B).sup.d] represents d sequential applications of the backward difference operator 1-B such that the original time series {[z.sub.t]} is self-differenced d times before application of the AR operator [phi](B).

As a shorthand, an ARIMA model that consists of an AR operator of order p, d backward difference operations, and a MA operator of order q is abbreviated as (p,d,q). A (p,0,0) model is referred to as an AR model of order p, or AR(p) in shorthand notation, and a (0,0,q) model is referred to as a MA model of order q (MA(q) in shorthand). A (p,0,q) model is referred to as an autoregressive-moving average model (ARMA(p,q) in shorthand) and a (0,d,q) model is referred to as an integrated moving average model (IMA(d,q) in shorthand). Finally, a random walk model is (0,1,0), while a random-walk-plus-uncorrelated-observation noise (RWPUN) model is (0,1,1).

In general, ARIMA models represent nonstationary time series. However, the time series that results from applying the backward difference operator d times to a (p,d,q) ARIMA process is a stationary ARMA(p,q) process. Typical constraints imposed on ARIMA models are that, when B is regarded as a complex variable, the polynomial in B representing the generalized autoregressive operator has zeros on or outside the unit circle ([absolute value of B]=1), whereas the polynomials representing the ordinary autoregressive operator and the moving average operator have zeros outside the unit circle. Additionally, the white noise variance, [[sigma].sup.2.sub.a], must be finite. These constraints ensure that a model is invertible and stationary (or, if not stationary, that at least it is of a suitable form).

One final property of ARIMA models is of interest here. The power spectrum for an ARMA (p,q) process represented by equation 1 is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (3)

where B = [e.sup.i2[pi]f; 0 [less than or equal to] f [less than or equal to] 1/2,

and the forward shift operator, F = [B.sup-1], has the opposite effect to B (i.e., [F.sup.m] [z.sub.t] = [z.sub.t+m]).

ARIMA models of time series with observation noise

Suppose, then, that [y.sub.1], [y.sub.2], ..., [y.sub.T] represent a time series of observations at times t = 1,2,... T such that

[y.sub.t] = [z.sub.t] + [e.sub.t], (4)

where the [z.sub.t]'s represent the unobserved, underlying process and the [e.sub.t]'s are IID normal variables with variance [[sigma].sub.e.sup.2] (i.e., e~N(0, [[sigma].sub.e.sup.2])) that are also independent of the [z.sub.t]'s (i.e., < [e.sub.j] [z.sub.k]>=0 for all j,k). The goal is to estimate the unobserved time series {[z.sub.t]} by using the observed time series {[y.sub.t]}. For the analysis of fishery-independent time series, it seems reasonable to assume that only the ARIMA model for the observed time series {[y.sub.t]} is known (it can be estimated using standard techniques). In particular, this assumption means that the model for {[z.sub.t]} is unknown within constraints implied by the observation equation. However, to develop the approach used here it is helpful to start as though the model for the unobserved process {[z.sub.t]} were known.

Thus, we assume that the time series {[z.sub.t]} can be represented by an ARIMA (p,d,q) process:

[phi](B)[z.sub.t] = [alpha](B)[c.sub.t], (5)

where the [c.sub.t]'s are IID and [c.sub.t] ~ N(0, [[sigma].sup.2.sub.c]). Substituting in [y.sub.t]-[e.sub.t] for [z.sub.t] and rearranging, one obtains

[phi](B)[y.sub.t] = [alpha](B)[c.sub.t] + [phi](B)[e.sub.t], (6)

which can be expressed as an ARIMA model for {[y.sub.t]} of the form

[phi](B)[y.sub.t] = [eta](B)[d.sub.t], (7)

where the [d.sub.t]'s are IID, [d.sub.t]~N(0, [[sigma].sup.2.sub.d]) and the MA operator is [eta](B). Thus, the generalized AR operator for the observed {[y.sub.t]} is identical to that for the unobserved {[z.sub.t]}. Furthermore, because [alpha](B) has order q and [phi](B) has order p+d, [eta](B) must be of order max(p+d,q). This requirement for [eta](B) constrains the order of potential ARIMA models that could describe the observed process: if p+d[less than or equal to]q then the observed process is also a (p,d,q) process, otherwise it is a (p,d,p+d) process.

In the more realistic situation where {[y.sub.t]} is observed, we can determine the generalized AR operator [phi](B) and MA operator [eta](B) for the observed process. Its ARIMA model will be order (P,D,Q), say, where the minimum possible value for Q is P+D. The model for the corresponding unobserved, underlying process {[z.sub.t]} will have order (p=P, d=D, q[less than or equal to]max(P+D,Q)) and its associated generalized AR operator will also be [phi](B). Furthermore, recognizing that [c.sub.t] and [e.sub.t] are independently distributed, it can be shown that the moving average operator for the unobserved process, [alpha](B), is additionally constrained by the ARIMA model for the observed process such that the following relationship must be satisfied:

[[sigma].sup.2.sub.c][alpha](B)[alpha](F) = [[sigma].sup.2.sub.d][eta](B)[eta](F) - [[sigma].sup.2.sub.e][phi](B)[phi](F). (8)

In this equation, [[sigma].sub.e.sup.2], [eta](B),and [phi](B) are known from the ARIMA model for the observed process, whereas [[sigma].sub.e.sup.2], [[sigma].sub.e.sup.2], and [alpha](B) are unknown.

In general, many combinations of [[sigma].sub.e.sup.2], [[sigma].sub.e.sup.2], and [alpha](B) will satisfy the equality. Defining an "acceptable model" for the unobserved process as one that, given the model for the observed process, [alpha]B) satisfies the previous equation and its zeros are on or outside the unit circle, Box et al. (1978) show that 1) for every given model of an observed process, at least one acceptable model for the unobserved process exists; 2) for a given model of an observed process, the possible values of [[sigma].sub.e.sup.2] are bounded; and 3) for a given model, every [[sigma].sub.e.sup.2] between 0 and the upper bound ([K.sup.*], say) determines a unique acceptable model. The upper bound on the observation error variance, [K.sup.*], is determined from the constraint that, for a model of the unobserved process to be acceptable, [[sigma].sub.c.sup.2][alpha](B)[alpha](F)[greater than or equal to] 0 everywhere on the unit circle (i.e., the power spectrum of the corresponding MA process is non-negative definite). Then, from equation 8, [K.sup.*] is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (9)

and is completely determined by the ARIMA model for the observed process. When [[sigma].sub.e.sup.2] = [K.sup.*], the variance of the added white noise is maximal, as will be the smoothing of the observed time series.

It is instructive to interpret Equations 8 and 9 in terms of constraints on the power spectra of {[z.sub.t]}, {[y.sub.t]}, and {[e.sub.t]}, although this interpretation is strictly correct only when {[z.sub.t]} and {[y.sub.t]} are stationary. Let [p.sub.z](f), [p.sub.y](f), [p.sub.e](f) denote the power spectra for {[z.sub.t]}, {[y.sub.t]}, and {[e.sub.t]}, respectively. Recalling the definition of the power spectrum (Eq. 3), Equation 8 on the unit circle can be easily recast (multiply both sides by 2/{[phi](B) [phi](F)}) as

[p.sub.z](f) = [p.sub.y](f) - [p.sub.e](f) (10)

because the power spectrum for white noise is constant with frequency: [p.sub.e](F)=2[[sigma].sup.2.sub.e]. Because power spectra are nonnegative definite and [p.sub.e](f) does not depend on f, the maximum possible observation noise variance [K.sup.*] corresponds to the minimum of [p.sub.y](f) over f (see Fig. 1). Thus, Equation 9 can be recast as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (11)

[FIGURE 1 OMITTED]

Note also that Equation 10 can be recast again as

[p.sub.z](f) = (1 - [p.sub.e](f) / [p.sub.y](f))[p.sub.y](f) = [[omega].sup.*](f)[p.sub.y](f), (12)

where [[omega].sup.*](f) represents the term in parentheses. It will be seen below that [[omega].sup.*](f) is identical to the polynomial of smoothing weights (Eq. 14) on the unit circle.

Returning to the estimation problem now, knowledge of [[sigma].sup.2.sub.e], then, together with the ARIMA model for the observed process, is sufficient to estimate {[z.sub.t]} (Box et al., 1978). When t is not close to the endpoints of the time series, the smoothed (maximum likelihood) estimate of [z.sub.t], which we denote as E([z.sub.t]|y), is a symmetric moving average filter of {[y.sub.t]} (Cleveland and Tiao, 1976; Box et al., 1978):

E([z.sub.t]|y) = [omega](B)[y.sub.t], (13)

such that

[omega](B = [[sigma].sup.2.sub.a]/[[sigma].sup.2.sub.d] [alpha](B)[alpha](F)/[eta](B)[eta](F) = [1 - [[sigma].sup.2.sub.e]/[[sigma].sup.2.sub.d] [phi](B)[phi](F)/[eta](B)[eta](F)], (14)

where the second equality follows from equation 8. Note that [omega](B) is identical to [[omega].sup.*](f) (Eq. 12) on the unit circle (B=[e.sup.-i2[pi]f)]), so that Equations 13 and 14 are also interpretable in terms of relations among power spectra. Equation 13 is equivalent to Equation 2 in Pennington (1985). It turns out (Cleveland and Tiao, 1976) that equation 14 is also valid for t near the endpoints of the time series; thus, one merely uses the ARIMA model for the observed time series to hindcast or forecast additional "observations" as needed. Box et al. (1978) describe a method for calculating the coefficients of w(B); because this reference may be difficult to obtain, we repeat their description in the Appendix.

In the case where the model for the underlying process is (0,1,0) (i.e., a random walk model), the model for the observed process is (0,1,1) (i.e., a RWPUN, model). For this case, Pennington (1985) noted that the value of the MA parameter [[eta].sub.1] of the observed process ([theta] in his notation) is

[[eta].sub.1] = [[sigma].sup.2.sub.e] / [[sigma].sup.2.sub.d], (15)

so that Equation 14 is completely determined by the ARIMA model for the observed process and the observation error variance, [[sigma].sup.2.sub.e], can be estimated.

Application of the noise reduction algorithm to bottom trawl survey data

We applied Box et al.'s (1978) noise reduction algorithm to 18 time series of abundance indices for finfish (nine speciesxtwo seasons; species are listed in Table 1) on Georges Bank in the northwest Atlantic. Time series for the fall survey spanned 40 years (1963-2002), and the spring time series spanned 36 years (1968-2003).

Stratified random bottom trawl surveys have been conducted annually on the northeastern continental shelf of the United States from Cape Hatteras to the Gulf of Maine in the fall since 1963 and in the spring since 1968 by the National Marine Fisheries Service (NMFS), Northeast Fisheries Science Center (NEFSC) (Azarovitz, 1981; Anonymous, 1988; Reid et al., 1999). We derived annual time series of estimated total abundance during the spring and fall seasons for nine finfish species from trawl survey catch records for the Georges Bank area. Reported catches (biomass) from strata encompassing Georges Bank were expanded on a per-tow basis by using species and length-specific corrections for catchability and coefficients from Edwards (1968) and Harley and Myers (2001). Annual stratified mean (expanded) catch-per-tow was then calculated for each species for both seasonal surveys. Finally, annual total population abundance was calculated by applying a swept area factor to the stratified mean catch-per-tow. Abundance indices corresponding to fall surveys spanned 1963-2002, and those corresponding to spring surveys spanned 1968-2003. Following Pennington (1985), we assumed that changes in catchability induced a lognormal error structure on the observed time series. However, because some time series (notably those for the schooling pelagic fish herring and mackerel) included zeroes, a ln(x) transformation was not applicable to all series. Thus, before further analysis, all time series were ln(x+1) transformed.

For each resulting time series, we used SAS (vers. 8, SAS Institute, Cary, NC) to identify candidate ARIMA models, to estimate parameters for these models, and to extend the time series by using forecasts and hindcasts before applying the smoothing algorithm. Candidate model structures for each time series were based on examination of empirical autocorrelation, partial autocorrelation, and inverse autocorrelation functions for the series (Box and Jenkins, 1976). When these functions indicated that the series was nonstationary, we applied the backward difference operator (1-B) to the series and examined the correlation functions for the new (differenced) series. In addition, because of its special significance in terms of interpretation, we always tried a RWUPN model as a candidate. Once candidate models were identified, we estimated coefficients for each model and calculated the associated Akaike information criterion (AIC; Akaike, 1973). AIC provides an objective criterion based on information theory for selecting the "best" approximating model from among a group of good candidate models (i.e., the criterion selects the model with the minimum AIC).

We examined the residuals and the empirical autocorrelation, inverse autocorrelation, and partial autocorrelation functions of the residual time series for significant deviation from white noise. When significant deviation was indicated, we dropped the model from further consideration. We also dropped models with orders (P,D,Q) that were inconsistent with the assumption of additive (after log transformation) observation noise (Q[greater than or equal to]P+D). Following this screening procedure, we were left with a group of "good" alternative models. We then selected the ARIMA model with the smallest AIC from among the remaining candidates as the "best" model to smooth the data.

We developed a MATLAB (vers. 6.5, The Mathworks Inc., Natick, MA) program to perform the noise reduction for each time series based on an extended version of the series and its associated ARIMA model. To allow smoothed estimates to be calculated near the ends of each time series, we extended each time series before smoothing to 40 years before its start by hindcasting with the selected ARIMA model and to 40 years past its end by forecasting with the model. Using the MATLAB program, we then calculated [K.sup.*], the maximal smoothing weights [omega](B) corresponding to [[sigma].sub.e.sup.2] = [K.sup.*], and the smoothed estimates E([z.sub.t]|y), following the approach of Box et al. (1978) outlined previously.

Results

The abundance indices we derived from fishery-independent bottom trawl surveys for nine finfish species during two seasons displayed a wide variety of trends, as well as a high degree of apparent interannual variability. This variability may be associated with environmentally driven, high-frequency changes in catchability, but we regarded it here as observation noise (Figs. 2 and 3). For example, springtime winter skate (Leucoraja ocellata) biomass appeared to increase by a factor of six during the early 1980s from a lower (but highly variable) mean level of ~140,000 t, to which it subsequently returned in the early 1990s (Fig. 2A). Yellowtail flounder (Limanda ferruginea), in contrast, declined by a factor of 10 during the 1970s and early 1980s from a high at the beginning of the time series and then began to rebound in the latter 1980s (Fig. 2G). Most recently, yellowtail flounder appear close to regaining their earlier peak abundance.

[FIGURES 2-3 OMITTED]

The 18 ARIMA models corresponding to the ln(x+1)-transformed abundance indices formed a diverse set (Tables 2 and 3). Half the time series were found to be nonstationary; however, one application of backward differencing (d=1) sufficed to achieve stationarity in all nine instances. Interestingly, the ARIMA models for all nine of these time series were consistent with RWPUN models (i.e., [0,1,1] models). For all the models, autoregressive orders (p) ranged between 0 and 3, and moving average orders (q) ranged from 1 to 5. The spring and fall models for little skate (L. erinacea) had the highest AR order, and the spring model for silver hake (Merluccius bilinearis) had the highest MA order. Only the models for little skate, Atlantic herring (Clupea harengus), and yellowtail flounder exhibited the same ARIMA order for both the spring and fall. And although none of the models reflected pure AR processes (q=0, an impossibility given the observation noise assumption), two were found to be pure MA processes (winter flounder and Atlantic mackerel [Scomber scombrus], both in spring). All nine IMA processes were RWPUN processes (i.e., (0,1,1)).

The effect of ARIMA-based noise reduction on these 18 sets of indices was fairly variable in outcome (Figs. 2 and 3). Little or no smoothing occurred for little skate (spring, Fig. 2B), silver hake (spring, Fig. 2D), and haddock (Melanogrammus aeglefinus) (fall, Fig. 3F). Moderate smoothing occurred for winter skate (spring, Fig. 2A; fall, Fig. 3A), little skate (fall, Fig. 3B), Atlantic herring (fall, Fig. 3C), haddock (spring, Fig. 2F), yellowtail flounder (spring, Fig. 2G; fall, Fig. 3G), winter flounder (Pseudopleuronectes americanus) (fall, Fig. 3H), and Atlantic mackerel (spring, Fig. 2I). A high degree of smoothing occurred for Atlantic herring (spring, Fig. 2C), silver hake (fall, Fig. 3D), Atlantic cod (Gadus morhua) (spring, Fig. 2E; fall, Fig. 3E), winter flounder (spring, Fig. 2H), and Atlantic mackerel (fall, Fig. 3I).

Impressions obtained from visually comparing the original and smoothed time series were consistent with the fraction of innovation variance ([K.sup.*]/[[sigma].sub.d.sup.2], Tables 2 and 3) in the observed time series ascribed to observation error. Recall that we chose to perform the maximal amount of smoothing possible (i.e., taking [[sigma].sub.d.sup.2]=[K.sup.*] in Eq. 13). Our visual classifications fell fairly neatly into the following categories based on [K.sup.*]/[[sigma].sub.d.sup.2]: 1) a low degree of smoothing corresponded to [K.sup.*]/[[sigma].sub.d.sup.2]<0.25, 2) a moderate degree of smoothing corresponded to 0.25< [K.sup.*]/[[sigma].sub.d.sup.2]<0.60, and 3) a high degree of smoothing corresponded to 0.60[greater than or equal to][K.sup.*]/[[sigma].sub.d.sup.2].

In addition, it appeared that the degree of smoothing varied inversely with the order of the moving average component of the ARIMA model for the observed time series. When the order of the moving average component was greater than 2, the degree of overall smoothing was typically small. Substantially more smoothing was evident when the order of the moving average component was 1, because of the equivalence with an exponentially weighted smoother.

The nine time series found to be RWPUN processes provided a means to check whether our choice to perform the maximal amount of smoothing was reasonable. For such time series, the value of the moving average coefficient, [[eta].sub.1], is equal to the ratio of the observation noise variance, [[sigma].sub.e.sup.2], to the innovation variance [[sigma].sub.d.sup.2]. Consequently, [[sigma].sup.2.sub.d]x[[eta].sub.1]=[[sigma].sup.2.sub.e]. For eight of the nine series, the ratio of [[sigma].sub.d.sup.2] to [K.sup.*] was greater than 80% (Table 4).

Discussion

Abundance indices derived from large-scale fishery-independent surveys typically exhibit interannual variability much higher than one would expect from within-survey variance (Pennington, 1985). Although true variability in the underlying population due to population-dynamic processes is reflected in the variability of an index, so too is observational noise arising both from within-survey sampling variability as well as from environmentally driven factors that affect catchability. Low signal-to-noise ratios in abundance indices due to high observational noise reduce one's ability to detect important changes or trends in actual population abundance.

To reduce the impact of white observation noise on time series data, Cleveland and Tiao (1976) developed an approach to noise reduction and smoothing that was based on their knowledge of the ARIMA model (Box and Jenkins, 1976) for an associated unobserved but underlying stochastic process. Box et al. (1978) extended this approach to address the situation where the ARIMA model for the underlying process was unknown, relying instead on an ARIMA model associated with the observed time series. In general, and certainly in regard to fishery-independent survey data, a model structure for the unobserved, underlying process will not be available. Hence, Box et al.'s (1978) approach will be the norm.

In the situation where the observed time series is stationary, we found that a frequency domain interpretation of Box et al.'s (1978) algorithm is particularly enlightening. When the times series is stationary, observation noise increases the power spectral density (PSD) of the observed process over that of the unobserved process by a fixed amount at all frequencies (Fig. 1). Consequently, the PSD of the unobserved process has the same shape as the PSD of the observed process, but with a fixed amount removed at all frequencies. Because the unobserved PSD must be nonnegative, the maximum observation noise (and thus the maximum noise reduction and smoothing that may occur) consistent with additive white noise corresponds to the minimum value of the observed process's PSD. Presumed higher values for the observation noise would require that the PSD for the unobserved process be negative over some frequency range. Box et al.'s (1978) algorithm computes the maximum possible observation noise and uses a time domain formula to calculate a smoothed, "unobserved" time series consistent with additive white noise for any noise level up to the maximum.

The ARIMA-based noise reduction approach was first applied to fisheries trawl survey data by Pennington (1985), who developed an alternative algorithm to that of Box et al. (1978) for estimating the smoothed time series. This algorithm was subsequently used in several studies (Pennington, 1985; Fogarty et al. (1); Pennington, 1986; Anonymous, 1988, 1993) to smooth time series of abundance indices derived from trawl survey data for the northeast coast of the United States. Of these studies, perhaps only Pennington (1986) constitutes a convincing demonstration of the utility of the ARIMA-based approach to time series noise reduction when the ARIMA model for the underlying process is unknown. Pennington's (1985) demonstration relied on an ARIMA model developed for the (usually unobserved) underlying population that was based on stock assessment results for the particular species considered. In contrast, RWPUN models were assumed a priori in Fogarty et al. (1) and Anonymous (1988, 1993) because the short time series (<25 observations per series) considered made identification of the underlying ARIMA models problematic. Only in Pennington (1986) was a model used that was fitted to trawl survey data and tested for appropriateness.

However, substantially longer time series are now available to test the ARIMA-based noise reduction concept. Consequently, to revisit the utility of the ARIMA-based approach for smoothing time series data derived from fisheries-independent trawl survey data, we applied Box et al.'s (1978) approach to smoothing to time series of annual abundance indices derived from the NMFS/NEFSC fisheries-independent bottom trawl survey of nine finfish species during two seasons on Georges Bank. Time series for the fall spanned 40 years (1963-2002), and the spring time series spanned 36 years (1968-2003). The species, comprising two elasmobranchs, three groundfish, two flatfish, and two pelagic schooling species, presented a broad range of life history characteristics (Table 1).

The noise reduction results we obtained varied among species and between seasons. Despite smoothing at the maximum level of noise reduction consistent with each model, very little smoothing was obtained for haddock (fall), little skate (spring), and silver hake (spring). The models for these time series had among the highest moving average orders (3-5). Examination of the PSD for each model in the frequency domain revealed very small minima, indicating no scope for noise reduction. Typically, models that had a MA order >2 exhibited substantially less smoothing than models with a MA order [less than or equal to]2. Models that were of MA order 1 generally resulted in the greatest smoothing. Models that were of MA order 0 did not (and could not) occur.

Of the 18 time series we considered (Tables 2 and 3, Figs. 2 and 3), only half were adequately represented as random-walk-plus-uncorrelated-noise (RWPUN) models. The ARIMA models we developed were varied in structure, ranging from a simple MA(1) model to rather complicated models with multiple parameters. Thus, our results provide evidence against the appropriateness of assuming a particular model structure a priori when the objective of the analysis is to identify the underlying dynamic structure of the population. This evidence is further strengthened by the results of Becerra-Munoz et al. (1999), who found only 9 of 52 abundance time series for finfish species from the NMFS/NEFSC bottom trawl survey that corresponded to random walk models.

As an exercise, we also attempted to smooth the nine data sets that were not adequately described as RWPUN models, using this model structure as an a priori assumption, even though our analysis indicated that other models were more appropriate. We were not able to estimate convergent models for three species: Atlantic cod (fall), winter flounder (spring), and Atlantic mackerel (spring). For the remaining six time series, the smoothed results appeared to be quite reasonable (Fig. 4), although we obtained little noise reduction when we employed the "correct" ARIMA model. The RWPUN-smoothed time series for haddock (Fig. 4, C and F) were similar to that for spawning biomass derived from virtual population analysis (see Brodziak et al. (2)), but the smoothed time series for silver hake (Fig. 4E) exhibited higher frequency variability than that found for total biomass with a production model (see Brodziak et al. (3)). From the standpoint of estimating the unobserved underlying process, these smoothed results should be viewed with some skepticism: the use of the RWPUN model is rather arbitrary in this situation and it may impose artificial structure on the smoothed results. However, it may be that these time series do not meet one of the key assumptions of the noise reduction method: namely that the observation noise is uncorrelated. The ARIMA models for all six time series had MA orders [greater than or equal to]3, and one effect of correlated observation noise could be to increase the MA order of the observed time series beyond that expected for uncorrelated noise.

[FIGURE 4 OMITTED]

Of the nine species considered, only three (little skate, Atlantic herring, and flounder) had models that exhibited the same ARIMA order for both the fall and spring surveys. Taken at face value, this would indicate that the other six species exhibited substantially different dynamical processes during the fall and spring that influenced abundance on Georges Bank. One potential mechanism for this could be differential seasonal migration patterns that result in changes in catchability that have different autocorrelation structures. For example, cod are distributed across the bank during the spring survey, but are found only in deeper waters on the periphery of the bank during the fall where they may be less available to the survey. Assuming that all Georges Bank cod are available to the spring survey, if the fraction of cod available to the survey in the fall is density dependent or is driven by autocorrelated environmental conditions, then the fall survey abundance will exhibit dynamical behavior different from that of the spring survey abundance (note: if this were the case, it would be inappropriate to smooth the fall time series for cod with the ARIMA noise reduction approach applied in our study because, as noted previously, one of the basic assumptions with this approach is that the observation noise is uncorrelated). Other plausible mechanisms can be developed, as well. However, we feel it more likely that the inconsistency in ARIMA order between spring and fall surveys for the same species is an estimation problem and indicates that even a 40-year time series may not be long enough to reduce the variability inherent in ARIMA model estimation to reasonable levels for most species.

Alternative methods, such as locally weighted scatterplot smoothing (LOESS), moving average filters, exponential smoothing filters, Kalman filters, and frequency-domain approaches can be applied to time series to achieve smoother results (e.g., Cleveland and Grosse, 1991; Hamilton, 1994). These approaches typically employ at least one user-determined parameter that can be used to change the amount of smoothing that an algorithm achieves. Generally, one "fiddles" with the adjustable parameters until a "nice," smoothed fit is achieved. However, we think it important to distinguish between these smoothing algorithms and the ARIMA-based noise reduction algorithms. It is quite possible to smooth out real fluctuations in the underlying population process. The principal advantage that we see for the ARIMA-based noise reduction algorithms (used with an appropriate model) over alternative methods is that the former provide a more objective approach to determining an appropriate level of smoothing. As noted previously, Pennington (1985) showed that when a RWPUN model is appropriate, the ARIMA smoothing approach is completely determined by the ARIMA model for the observed time series because it is possible to determine the observation noise variance from the model parameters. In the more general case, Box et al.'s (1978) algorithm at least yields a maximum value for the variance of the observational noise and thus sets an upper limit to the amount of noise reduction and smoothing that can be achieved. For trawl survey data, our results from nine time series where RWPUN models were appropriate (and we can consequently estimate the actual observation noise variance) indicate that smoothing at ~90% of the maximum possible noise reduction level is not an unreasonable default percentage (Table 4).

One drawback to the greater application of ARIMA-based noise reduction methods to time series data is the lack of an integrated software package that allows a user 1) to quickly evaluate an appropriate ARIMA model for a given time series, and 2) to calculate the smoothed time series. We used SAS for the first step and MATLAB for the second, but we found this arrangement rather awkward and burdensome. However, econometrically oriented software packages such as ForecastPro (4) or AutoBox (5) that automate model selection may substantially simplify the first step even if they don't address the second step.

On the whole, ARIMA-based time series models appear to provide the basis for a more objective approach to reducing observation noise in time series data, including time series of fishery abundance indices derived from trawl survey data, than do more conventional smoothing approaches. In the absence of additional information regarding the level of observation noise, we recommend smoothing trawl survey data at 90% of the maximum possible noise reduction level. We also suggest that development of an integrated software package for implementing ARIMA-based noise reduction will facilitate future use of this method.

Finally, if a smoothed time series is desired (e.g., for graphical presentation only), then use of a RWPUN model in lieu of a model-fitting exercise will generally yield a curve pleasing to the eye. Alternatively, other methods such as LOESS could be employed to generate the smoothed results. However, if the resulting time series is to be used for further analysis of the dynamical behavior of the fish stock, we strongly recommend that a model-fitting approach be used to identify the most appropriate ARIMA model for the observed time series, from which the time series for the unobserved, underlying process can be computed. Otherwise, real fluctuations in the underlying process may be oversmoothed, resulting in an apparent dynamical behavior that displays little variability. This oversmoothing, in turn, may lead to erroneous conclusions being drawn regarding, for example, the resiliency of a stock to exploitation or environmental change, and to perhaps concomitant errors being propagated in advice provided to fishery managers.

Appendix

For the convenience of the reader, we summarize here Box et al.'s (1978) algorithm to calculate the coefficients of the smoothing polynomial [omega](B). Recall from Equation 14 that

[omega](B) = [1 - [[sigma].sup.2.sub.e]/[[sigma].sup.2.sub.d] [phi](B)[phi](F)/[eta](B)[eta](F], (A.1)

where [phi](B) is a polynomial of order P+D and [eta](B) is a polynomial of order Q. Because Q [greater than or equal to] P+D, one can write

[phi](B) = 1 - [[phi].sub.1][B.sup.1] - [[phi].sub.2][B.sup.2] ... - [[phi].sub.Q][B.sup.Q]; [eta](B) = 1 - [[eta].sub.1][B.sup.1] - [[eta].sub.2][B.sup.2] ... - [[eta].sub.Q][B.sup.Q], (A.2)

where [[phi].sub.j] = 0 for j>P+D.

First, define

C(B) [equivalent to] [phi](B)/[eta](B); C(B) = 1 + [C.sub.1][B.sup.1] + [C.sub.2][B.sup.2] +.... (A.3)

One can solve for the coefficients of C using an iterative process by recognizing that the coefficients of each power of B in the following expression must be zero:

0 = [phi](B)- [eta](B)C(B).

Consequently, one obtains

[C.sub.1] = [[eta].sub.1] - [[phi].sub.1]; [C.sub.j] = [[eta].sub.j] - [[phi].sub.j] + [j-1.summation over (i=1)] [C.sub.j-i][[eta].sub.i], j = 2,... Q (A.4)

and higher order coefficients of C (i.e., for j>Q) can be computed recursively using the relation

[C.sub.j] = [Q.summation over (i=1)] [C.sub.j-i][[eta].sub.i]. (A.5)

Next, define

X(B,F) [equivalent to]C(B)[[phi](F)/[eta](F) (A.6)

where X(B,F) = [X.sub.0] + [[infinity].summation over (j=1)] [X.sub.j] ([B.sup.j] + [F.sup.j]).

From this definition, one obtains

X(B, F)[eta](F) [equivalent to]C(B)[phi](F). (A.7)

The largest degree of F on the righthand side of this equation is Q. By matching coefficients of [F.sup.j] on both sides of the equation, Box et al. (1978) found, for j=0,1,... Q, that

[eta]X = C[phi], (A.8)

where X, [phi] are Q+1 column vectors [X.sup.T]=([X.sub.0], [X.sub.1],..., [X.sub.Q]), [[phi].sup.T]=(1, -[[phi].sub.1],..., -[[phi].sub.Q]) and [eta], C are (Q+1)x(Q+1) matrices such that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A.9)

Thus, Equations A.8 and A.9 define a set of Q+1 equations in Q+1 unknowns ([x.sub.0], [X.sub.1],..., [X.sub.Q]) which can be solved for by matrix inversion such that

X = [n.sup.-1]C[phi]. (A.10)

For j>Q, the [X.sub.j]'s are computed recursively by using the relation

[X.sub.j] = [Q.summation over (i=1)] [X.sub.j-1][[eta].sub.i]. (A.11)

Finally, the coefficients of [omega](B) are given by

[[omega].sub.0] = 1 - [[sigma].sup.2.sub.e]/[[sigma].sup.2.sub.d] [X.sub.0]; [[omega].sub.j] = -[[sigma].sup.2.sub.e]/[[sigma].sup.2.sub.d] [X.sub.j], j = 1,.... (A.12)

Acknowledgments

We would like to thank M. Pennington for his insightful comments on an early version of this manuscript. The comments and suggestions from two anonymous reviewers were also very helpful and are much appreciated. This work was supported by the National Research Council, which provided funding support for one of us (Stockhausen) as a postdoctoral fellow at the Northeast Fisheries Science Center.

Manuscript submitted 12 November 2004 to the Scientific Editor's Office.

Manuscript approved for publication 18 April 2006 by the Scientific Editor.

Literature cited

Akaike, H. 1973. Information theory as an extension of the maximum likelihood principal. In Second international symposium on information theory (B. N. Petrov, and F. Csaki, eds.), 267-281. Akademiai Kiado, Budapest, Hungary.

Anonymous. 1988. An evaluation of the bottom trawl survey program of the Northeast Fisheries Center. NOAA Tech. Memo. NMFS-F/NEC-52, 83 p. Northeast Fish. Center, 166 Water St., Woods Hole, MA 02543.

1993. Status of fishery resources off the northeastern United States for 1993. NOAA Tech. Memo. NMFSF/NEC-101, 140 p. Northeast Fish. Center, 166 Water St., Woods Hole, MA 02543.

Azarovitz, T. R. 1981. A brief historical review of the Woods Hole Laboratory trawl survey time series. In Bottom trawl surveys (W. G. Doubleday, and D. Rivard, eds.), 62-67. Can. Spec. Pub. Fish. Aquat. Sci., vol. 58.

Becerra-Munoz, S., D. B. Hayes, and W. W. Taylor. 1999. Stationarity and rate of damping of modeled indices of fish abundance in relation to their exploitation status in the Northwest Atlantic Ocean. Ecol. Mod. 117:225-238.

Box, G. E. P., S. C. Hillmer, and G. C. Tiao. 1978. Analysis of seasonal time series. In Seasonal analysis of economic time series (A. Zellner, ed.), 309-334. U.S. Dep. Comm., Washington, DC.

Box, G. E. P., and G. M. Jenkins. 1976. Time series analysis: forecasting and control, 575 p. Holden-Day, Oakland, CA.

Byrne, C. J., T. R. Azarovitz, and M. P. Sissenwine. 1981. Factors affecting variability of research trawl surveys. Can. Spec. Pub. Fish. Aquat. Sci. 58:238-273.

Clark, S. H. 1979. Application of bottom trawl survey data to fish stock assessment. Fisheries 4:9-15.

Cleveland, W. S., and E. Grosse. 1991. Computational methods for local regression. Stats. and Comp. 1:47-62.

Cleveland, W. P., and G. C. Tiao. 1976. Decomposition of season time series: a model for the Census X-11 program. J. Am. Stat. Assoc. 71:581-587.

Collie, J. S., and M. P. Sissenwine. 1983. Estimating population size from relative abundance data measured with error. Can. J. Fish. Aquat. Sci. 40:1971-1983.

Edwards, R. L. 1968. Fishery resources of the North Atlantic area. In The future of the fishing industry in the United States (D. W. Gilbert, ed.), p. 52-60. Univ. Washington, Seattle, WA.

Enders, W. 2004. Applied econometric time series, 460 p. John Wiley and Sons, Inc., Hoboken, NJ.

Grosslein, M. D. 1969. Groundfish survey program of BCF Woods Hole. Comm. Fish. Rev. 31:22-30.

Hamilton, J. 1994. Time series analysis, 820 p. Princeton Univ. Press, Princeton, NJ.

Harley, S. H., and R. A. Myers. 2001. Hierarchical Bayesian models of length-specific catchability of research trawl surveys. Can. J. Fish. Aquat. Sci. 58:1569-1584.

Helser, T. E., and D. B. Hayes. 1995. Providing quantitative management advice from stock abundance indices based on research surveys. Fish. Bull. 93:290-298.

Lloret, J., J. Lleonart, and I. Sole. 2000. Time series modeling of landings in northwest Mediterranean Sea. ICES J. Mar. Sci. 57:171-184.

Nogueira, E., F. F. Perez, and A. F. Rios. 1998. Modelling nutrients and chlorophyll a time series in an estuarine upwelling ecosystem (Ria de Vigo: NW Spain) using the Box-Jenkins approach. Estuar. Coast. Shelf Sci. 46:267-286.

Pennington, M. 1985. Estimating the relative abundance of fish from a series of trawl surveys. Biometrics 41:197-202.

1986. Some statistical techniques for estimating abundance indices from trawl surveys. Fish. Bull. 84: 519-525.

Reid, R. N., F. P. Almeida, and C. A. Zetlin. 1999. Fishery-independent surveys, data sources, and methods. NOAA Tech. Memo. NMFS-NE-122. Northeast Fish. Sci. Center, 166 Water St., Woods Hole, MA 02543.

(1) Fogarty, M. J., J. S. Idoine, F. P. Almeida, and M. Pennington. 1986. Modeling trends in abundance based on research vessel surveys. ICES CM (council meeting) 1986/G, p. 92. ICES, Copenhagen, Denmark.

(2) Brodziak, J., M. Traver and L. Col. 2005. Georges Bank haddock. In Assessment of 19 northeast groundfish stocks through 2004 (R. K. Mayo, and M. Terceiro, eds.), section 2, p. 30-80. 2005 groundfish assessment review meeting. Northeast Fisheries Science Center, Woods Hole, Massachusetts; 15-19 August 2005. NEFSC Ref. Doc. 05-13. NEFSC, 166 Water Street, Woods Hole, MA 02543.

(3) Brodziak, J. K. T., E. M. Holmes, K. A. Sosebee, and R. K. Mayo. 2001. Assessment of the silver hake resource in the Northwest Atlantic in 2000, 134 p. NEFSC Ref. Doc. 0103. NEFSC, 166 Water Street, Woods Hole, MA 02543.

(4) Business Forecast Systems, Inc. 2006. Website: http://www. forecastpro.com/products/fpfamily/index.html (accessed on 29 March 2006).

(5) Automatic Forecasting Systems. 2003. Website: http:// www.autobox.com/autoboxdesc.htm (accessed on 29 March 29 2006).

William T. Stockhausen (contact author)

Michael J. Fogarty

Email address for W. T. Stockhausen: William.Stockhausen@noaa.gov

National Ocean and Atmospheric Administration

National Marine Fisheries Service

Northeast Fisheries Science Center

166 Water Street

Woods Hole, MA 02543

Present address for corresponding author: National Marine Fisheries Service

Alaska Fisheries Science Center

7600 Sand Point Way NE

Seattle, Washington 98115

Table 1 Time series of abundance indices for the following finfish species on Georges Bank were derived from a fisheries-independent trawl survey and used to test the ARIMA-based smoothing algorithm. Common name Scientific name Type Winter skate Leucoraja ocellata Elasmobranch Little skate Leucoraja erinacea Elasmobranch Silver hake Merluccius bilinearis Groundfish Atlantic cod Gadus morhua Groundfish Haddock Melanogrammus aeglefinus Groundfish Winter flounder Pseudopleurenectes americanus Flatfish Yellowtail flounder Limanda ferruginea Flatfish Atlantic herring Clupea harengus Pelagic schooling fish Atlantic mackerel Scomber scombrus Pelagic schooling fish Table 2 ARIMA smoothing results for ln(x+l) transformed time series based on spring surveys. Here [[sigma].sub.d.sup.2] ad is the innovation variance for the transformed time series and [[kappa].sup.*] is the maximum consistent observation noise variance, scaled to [[sigma].sub.d.sup.2] (i.e., [[KAPPA].sup.*] [[sigma].sub.d.sup.2]). The degree of smoothing is based on a visual comparison of the original and smoothed time series. Species ARIMA model (p,d,q) [[sigma].sub.d.sup.2] Winter skate (0, 1, 1) 0.33 Little skate (3, 0, 4) 0.30 Atlantic herring (0, 1, 1) 4.90 Silver hake (0, 0, 5) 0.85 Atlantic cod (0, 1, 1) 0.47 Haddock (1, 0, 3) 0.42 Yellowtail flounder (0, 1, 1) 0.45 Winter flounder (0, 0, 1) 0.33 Atlantic mackerel (0, 0, 2) 16.9 Species [[kappa].sup.*] Degree of smoothing Winter skate 0.52 medium Little skate 0.14 low Atlantic herring 0.59 high Silver hake 0.02 low Atlantic cod 0.71 high Haddock 0.29 medium Yellowtail flounder 0.49 medium Winter flounder 0.89 high Atlantic mackerel 0.34 medium Table 3 ARIMA smoothing results for ln(x+1) transformed time series based on fall surveys. Here, [[sigma].sub.d.sup.2] is the innovation variance for the transformed time series and [[kappa].sup.*] is the maximum consistent observation noise variance, scaled to [[sigma].sub.d.sup.2] (i.e., [[KAPPA].sup.*] [[sigma].sub.d.sup.2]). The degree of smoothing is based on a visual comparison of the original and smoothed time series. Species ARIMA model (p,d,q) [[sigma].sub.d.sup.2] Winter skate (2, 0, 3) 0.17 Little skate (3, 0, 4) 0.12 Atlantic herring (0, 1, 1) 4.75 Silver hake (0, 1, 1) 0.22 Atlantic cod (1, 0, 1) 0.55 Haddock (2, 0, 3) 0.75 Yellowtail flounder (0, 1, 1) 0.40 Winter flounder (0, 1, 1) 0.40 Atlantic mackerel (0, 1, 1) 7.00 Species [[kappa].sup.*] Degree of smoothing Winter skate 0.31 medium Little skate 0.38 medium Atlantic herring 0.48 medium Silver hake 0.61 high Atlantic cod 0.77 high Haddock <0.01 low Yellowtail flounder 0.45 medium Winter flounder 0.50 medium Atlantic mackerel 0.84 high Table 4 For time series consistent with random-walk-plus-uncorrelated-noise models, the table provides a comparison between the estimated observation noise variance ([[sigma].sub.e.sup.2]) and the maximum noise variance ([[KAPPA].sup.*]). The latter was used to smooth all time series. Season Species [n.sub.1]= [[sigma]. [[sigma].sup.2.sub.e] sup.2.sub.d] [[sigma].sup.2.sub.d] Fall Atlantic herring 0.392 4.75 Silver hake 0.566 0.22 Yellowtail flounder 0.342 0.40 Winter flounder 0.420 0.40 Atlantic mackerel 0.832 7.00 Spring Winter skate 0.438 0.33 Atlantic herring 0.541 4.90 Atlantic cod 0.680 0.47 Yellowtail flounder 0.399 0.45 Season Species [[KAPPA].sup.*] [[sigma]. sup.2.sub.e] [[KAPPA].sup.*] Fall Atlantic herring 2.28 0.82 Silver hake 0.13 0.93 Yellowtail flounder 0.18 0.76 Winter flounder 0.20 0.84 Atlantic mackerel 5.88 0.99 Spring Winter skate 0.17 0.84 Atlantic herring 2.89 0.92 Atlantic cod 0.33 0.96 Yellowtail flounder 0.22 0.81

Printer friendly Cite/link Email Feedback | |

Title Annotation: | autoregressive integrated moving average |
---|---|

Author: | Stockhausen, William T.; Fogarty, Michael J. |

Publication: | Fishery Bulletin |

Geographic Code: | 1USA |

Date: | Jan 1, 2007 |

Words: | 9381 |

Previous Article: | Ecotypic variation and predatory behavior among killer whales (Orcinus orca) off the eastern Aleutian Islands, Alaska. |

Next Article: | Abundance and distribution of the eastern North Pacific Steller sea lion (Eumetopias jubatus) population. |

Topics: |