Climate variability occurs across a spectrum of spatial (global, regional, local) and temporal (interdecadal, interannual, seasonal) scales. These interacting processes combine to form a complex picture that often defies simple interpretation. Marine fish populations respond to climate variability through a variety of physical and biological mechanisms. Traditional population dynamics models have generally ignored the effect of climate on fish population fluctuations, instead relating population size to biological (stock size, fecundity) and harvesting factors. However, in the hope of improving management practices, more research needs to be directed towards better understanding the relationship between environmental variability and fish production dynamics (Shepard et al. 1984).
The study of climate change and fish populations requires going beyond the traditional statistical tools employed in fisheries research. A few of the reasons for employing better statistical tools are:
1) Analyses of climate change impact fall under the historical/comparative scientific method (Gould 1989, Francis and Hare 1994) rather than the experimental/predictive. Without controlled and replicated situations the suite of confirmatory statistical methods (e.g., analysis of variance, t test) are not applicable.
2) In such analyses, one is examining variability over large spatial scales. A single point variable, such as sea surface temperature at Ocean Station P or sea level pressure at Dutch Harbor may not be representative of the larger scale variability, and may even be misleading. Thus, one is led to analyses that explicitly recognize the larger scale patterns of variability.
3) Since climate change analyses involve observations over time, one will likely encounter the effects of autocorrelation, or lack of independence between successive data points. The presence of autocorrelation greatly influences the significance of trends in individual series, as well as the appearance of a relationship between two time series.
In this chapter, I describe the motivation, methodology and notation for two statistical methodologies that address the concerns cited above, and are heavily used in the analyses conducted in this study. The first methodology, principal component analysis (PCA), has been widely used in climatological research to identify patterns of climate variability and how they evolve over time. PCA is closely related to the other "eigentechniques": empirical orthogonal function analysis and common factor analysis. I develop the PCA methodology in detail and note similarities and differences with respect to the other eigentechniques. I also summarize recent advances and controversies in PCA concerning the issues of rotation and eigenvector interpretability. These issues strongly affect the outcome of PCA and have received little, if any, attention in the fisheries literature. Secondly, I review the relevant aspects of time series analysis, including auto and cross correlation functions, prewhitening, Box-Jenkins AutoRegressive Integrated Moving Average (ARIMA) models, and Intervention Analysis. A major goal of this review is to call attention to fisheries and climate-specific issues that must be addressed in these type of analyses and provide an introduction to the sizable statistically-oriented climate related literature that has developed over the past two decades.
The goal of PCA, in a climatological context, is to redistribute the variability of a large set of variables to a much smaller set of variables (components) that 1) contain most of the original variance and 2) can be interpretable in terms of physically meaningful processes. For example, suppose one wished to analyze the historical variability of sea surface temperature (SST) in the North Pacific Ocean. Commonly available data sets of North Pacific SST can be quite large: in the example that follows, I analyze 129 months (43 years of three month winters) of data at 100 grid points. Each grid point has a unique temporal history, with varying degrees of similarity to other grid points. One might surmise that each point's history reflects its exposure to large-scale, possibly overlapping, physical processes (e.g., seasonal cycle, El Niños, winter storms, climatic regimes) plus some amount of measurement error and random noise. An economical means of describing SST variability would be to have estimates of the influence of these processes for each grid point. El Niño events, for instance, could be represented by a map showing the range (and relative regional strength) of influence plus a time series that documents the occurrence and relative temporal strength of each event. What PCA does, in summary, is to estimate both the physical process maps (spatial variability) and time series coefficients modulating each process (temporal variability) from the variance-covariance matrix of the observed data.
Situations such as that described above, and many far more complicated, frequently arise in climatological research. For that reason, PCA, and its relatives, empirical orthogonal function analysis (EOFA) and common factor analysis (CFA), have become standard analytical tools of meteorologists, oceanographers and climatologists. The use of eigentechniques in climatological research has been routine since the 1970s. Richman (1986) mentions that more than 80 PCA applications appeared in the meteorological literature between 1980 and 1986. Fukuoka (1951) appears to be the first to have applied PCA in a meteorological setting. However, it was Lorenz (1956) and Kutzbach (1967) who, by providing lucid explanations and demonstrations of the technique, were most responsible for the widespread adoption of PCA. Davis (1976) introduced the technique to the oceanographic community.
Before further description of PCA, it might be useful to clarify certain problems relating to eigentechnique terminology. Though often used interchangeably in the literature, PCA, EOFA, and CFA each represent a distinct "eigenmodel" whose specifications must be known in order to be interpretable. In their initial solutions, PCA and EOFA differ only by a scaling constant, which can become important when alternative solutions ("rotation" discussed below) are considered. Though mathematically similar to the other two techniques, CFA is a conceptually different type of analysis. Preisendorfer (1988) summarizes PCA as a tool for "isolating the sources of data variability", while CFA is characterized as appropriate for "studying sources of data covariability." Essentially, PCA and EOFA concentrate on the diagonal elements of the variance-covariance matrix, CFA concentrates on the off-diagonal elements (Jobson 1992). Throughout this exposition, and in the example, I will use the PCA model and terminology.
PCA input
The input to PCA consists of variables, cases and fixed entities. In a climatological setting, the candidates for analysis are station (a physical location such as latitude/longitude grid point or a city), time (usually in equidistant intervals) and field (climatological variable such as sea level pressure or geopotential height). A PCA is generally conducted on an "anomalous" (or "centered") dataset where the fixed entity mean has been removed from the case values of each variable. There are six basic ways (termed O-mode through T-mode) of specifying a PC analysis depending on what one chooses to be the variables, cases and fixed entities (Table 1, Figure 1.1). The large majority of climatological studies are S-mode analyses, with variables referencing stations, cases referencing time and fixed entity being the climatological field of interest. The example contained herein is an S-mode analysis in which I analyze how the sea surface temperature field varies (from the long-term mean field) over space (the North Pacific) and time (springs from 1950-1992). I will also maintain the S-mode terminology throughout my exposition of the technique. Examples of the other PC analysis modes are provided in Richman (1986) and Preisendorfer (1988).
PCA output
The output of PCA consists of three types of results: principal components (or "scores"), eigenvectors (also termed EOFs or "loadings") and eigenvalues. There are as many PCs, EOFs and eigenvalues as there are input original variables. EOFs and PCs are always "paired up", as in the linear model described above. In an S-mode analysis, the EOFs describe spatial variability, the PCs temporal variability. The EOFs are generally plotted as contour or vector maps, from which one can assess which regions are closely related, inversely related or unrelated, as well as identify centers of activity, regions with strong gradients, etc. The PCs, plotted as time series, quantify the overall strength of the associated EOF pattern over time. That is, the relative relationships between points on the grid (shown by the EOF) remain the same, but the absolute magnitude of the pattern changes with time (shown by the PC). The EOF/PC pairs (hereafter, "modes") are ranked by their importance to the overall variability of the dataset. The relative importance of each mode is determined by its associated eigenvalue, which is used to calculate the variance attributable to that mode. In a dataset consisting of one or more strong signals (or physical forcing processes), most of the original variability is captured by the first few PC/EOF modes.
In the initial PC solution, each EOF extracted following the first is orthogonal to all those preceding it. Successive PCs also carry the property of orthogonality. For this reason, PCA is often used to transform a large set of correlated variables to a smaller set of uncorrelated variables for use in other analyses such as multiple regression, while still retaining most of the original variance. The spatial orthogonality constraint can result in predictable geometric relations between the first EOF and subsequent ones, particularly for regularly gridded data (Horel 1981). Since one of the major goals of PCA is to recover spatial patterns that have a clear physical interpretation, it is often necessary to remove the orthogonality condition. Techniques for doing this are termed rotation, of which there are two types: orthogonal, which retains the orthogonality of the PCs while removing it from the EOFs, and oblique, which removes orthogonality from both. Much current literature (see Richman 1986 for a literature review) suggests that rotation is a necessary component of PCA (termed RPCA) and I discuss the various techniques and demonstrate an orthogonal solution, known as Varimax, in the section on Interpretation (section 2.2.4.2).
In its simplest form, PCA is used to detect the presence of "standing oscillations" in a single scalar field dataset. There have been a number of extensions to PCA. In combined PCA, one "chains" together fields of interest to isolate patterns of variability that occur synchronously (Kutzbach 1967). Thus, in addition to separate analyses of SLP and SST, one could combine them into a single PCA to determine if modes identified in separate analyses are temporally correlated. Also, combined PCA can be used to analyze vector fields such as surface winds (Barnett 1977). A form of combined PCA known as extended EOF analysis is used to study propagating features of climatological fields by chaining together at each time interval several time sequences of observations (Weare and Nasstrom 1982). A different form of PCA on vector valued fields, known as complex PCA, with the two components represented as the real and imaginary parts of complex numbers, is used to identify rotational variability in addition to spatial centers of action (Hardy 1977, Hardy and Walton 1978). Detection of propagating features can also be studied using the techniques of complex PCA (Barnett 1983, Horel 1984) or frequency domain PCA (Wallace and Dickinson 1972).
Additionally, several other techniques, not necessarily derived from PCA, have been applied to the study of single and multiple field variability. For examining coupled modes of variability between two fields, the techniques of canonical correlation analysis (Glahn 1968, Barnett 1981) and singular value decomposition (Prohaska 1976, Bretherton et al. 1992) have come into common use. The most recent advances combine EOF spatial representations with autoregressive moving average modeling of dynamic systems. These techniques are known as principal interaction patterns (PIP), principal oscillation patterns (POP), and complex POP analyses (Hasselmann 1988, Bürger 1993).
For such a vast subject, the basic mathematics of PCA are surprisingly straightforward. I will give just a summary of the relevant equations; for a full derivation see any of a number of texts on multivariate statistical methods (e.g., Jolliffe 1986, Jobson 1992). In essence one extracts the eigenvalues and eigenvectors of a square covariance or correlation matrix formed from the input data. To proceed, consider an S-mode analysis of a data field defined over N time steps on n grid points. The input data matrix X is of dimension N x n, with i indexing time (cases) and j indexing grid points (variables). The value of the field at any point in time and space is then denoted xij. In the basic PC model, a new set of variables Z are formed as linear combinations of the original variables, weighted by the coefficients E. It is hoped that most of the original variability of the n original variables will be concentrated in a much smaller number (r £ n, indexed by m) of transformed variables.
(1)
or, in matrix notation:
(2) 
The E matrix contains what will be subsequently referred to as the eigenvectors (or EOFs) and the Z matrix the PCs. In this arrangement, a single column of Z is one PC and a single row of E is its associated EOF. When rearranged as (3), we see that the original dataset can be reconstructed from the PCs and EOFs.
(3) 
The EOF coefficients are chosen so as to maximize each corresponding PC's variance, with the restriction that each PC be uncorrelated with previous PCs. The first principal component will account for the greatest amount of total variance, and each subsequent PC accounts for the maximum possible of the remaining variance. The sum of the squares of the EOF coefficients is constrained to unity ("unit length") since otherwise the PC variances could be made arbitrarily large: it is the relative magnitudes of the coefficients that is of interest. This also ensures that the total variance of the transformed and original datasets is equal. The EOFs are also uncorrelated with each other, hence the EOF matrix multiplied by its transpose yields the identity matrix:
(4) 
The variance/covariance matrix for the original (C) and transformed (L) datasets are given respectively by (5) and (6):
(5) 
(6) 
The following relationships then exist between C and L:
(7) 
(8) 
The last term in (8) derives from the fact that the off-diagonal elements of L are zero due to the orthogonality of the PCs. The solution to maximizing (2), i.e., the variance of successive PCs, subject to the condition (4) requires solving the matrix equation:
(9) [C - lI]E = 0
That is, one must find the characteristic roots (eigenvalues (l)) and eigenvectors of the data covariance matrix. There will be as many roots as original variables, however, many of these will not be significantly different from zero, leading to the reduction in the number of transformed variables referred to earlier. If the number of new variables formed, r, is equal to the original number of variables, n, the total variance of the dataset is equal to the total variance of the principal components, which is also the sum of the eigenvalues.
(10) 
and the amount of variance explained by each PC/EOF mode is calculated by dividing the corresponding eigenvalue by the trace of the eigenvalue matrix. In practice, the number of new variables, r, is usually much less than the original number, n. When a smaller set of transformed variables is selected, the fraction of total variance contained in the subset is given by the ratio of the traces of the covariance matrices. If a correlation dispersion matrix is used, the sum of the variance is equal to the number of variables (e.g., grid points).
Two techniques are generally used to extract the eigenvalues and eigenvectors of the covariance matrix C: the "power method" and the method of singular value decomposition (Jolliffe 1986). I use the SVD method in all of my examples and results. Both procedures are routinely available in computer software. A number of PC based packages include a PCA routine that includes many of the various options discussed below.
One final note concerning decomposition of the covariance matrix deserves mention. In an S-mode analysis, the covariance matrix is formed from the individual (time) values of the variables (grid points). If the number of variables (grid points) is larger than the number of cases (n > N), the resultant covariance matrix (n x n) will be even larger than the data matrix (N x n), perhaps resulting in a problem too costly in computing terms. One option is to form the covariance matrix after "swapping" the variables and cases. So, instead of forming the covariance matrix from the time values at each station, one forms the matrix from the station values at each time interval. This is equivalent to a T-mode analysis. The resultant covariance matrix is of size N x N, and the E matrix now represents the PCs and the Z matrix the EOFs. The same swapping of variables and cases can be done among other pairs of modes (Table 1.1) in order to work with the smaller covariance matrix. If, however, rotation is to be performed on the results based on a particular mode, the swapping cannot be done. For example, in an S-mode analysis, rotation is used to simplify the spatial interpretation of the field by grouping similar grid points. Using the T-mode covariance matrix, the effect of rotation would be to group similar time periods among the grid points.
In an S-mode analysis, the climatological mean is removed to generate t-centered data (i.e., time centered). At each grid point (Figure 1.2) a separate climatological mean is calculated. Note that there are different definitions of the climatological mean. For example, if consecutive monthly data over a number of years is analyzed, one has the option of removing from the monthly value at each grid point either 1) the corresponding long term monthly mean or 2) the long term annual mean across all 12 months. Most PC analyses are of type 1. The principal difference between the two is that the 2nd option includes an analysis of the annual cycle which may tend to dominate the analysis (Trenberth and Paolino 1981).
In the spring SST example, it was necessary to compute a seasonal average. To do so, I first computed annual spring SST values by averaging the monthly values for March, April, and May during each year. The long term mean spring SST field was than computed by averaging across the years 1950-1992. The seasonal mean field that was subtracted from each of the yearly seasonal maps of spring SST is illustrated in Figure 1.3.
Depending on the type of analysis, there are at least two situations where it may be necessary or advantageous to weight the input data. For example, in an S-mode analysis using a regular latitude/longitude grid there will be more stations per unit area in higher latitude regions (Figure 1.2, lower panel). If there is some level of spatial correlation (which is virtually always present), the points closer together will tend to show more similar patterns of variability thus biasing the covariance matrix and resultant eigenvectors. In such a situation, the "data rich" areas will dominate the results. There are two means of overcoming this problem. The first is to weight the input data by the square root of the cosine of the latitude of each station (North et al. 1982a). This has the effect of weighting the variance for each station by the size of the geographical region it represents. The relative size of the grid points in the bottom panel of Figure 1.2 illustrate the relative latitude weights. A map of the spring SST variance field (prior to weighting for PCA) is illustrated in Figure 1.3. Note that a reverse latitude weighting must later be applied to the transformed variables (eigenvectors) so that the correct scale is preserved. This is the method I used in all analyses. A second method, which also applies to irregularly spaced data, is to interpolate the data to an equal area grid (Karl et al. 1982).
As noted earlier, one can include multiple fields in a single PCA analysis. The fields are likely to have different units and levels of variability (e.g., SLP and SST). One solution is to use standardized values, however, this is equivalent to 1) using a correlation as opposed to covariance matrix (discussed below) and 2) partitioning the total variance equally between the two fields. An alternative is to weight the input variables by a predetermined fraction of its field's total variance to the unweighted data's total variance. For example, to equally weight two input fields without normalizing (and thereby losing real point to point differences in variance), one can multiply each input variable by the inverse of the ratio of the square root of its field's total variance to the square root of the input data total variance. Other weighting techniques have been suggested by Lorenz (1956).
The decision of whether to normalize input data and which dispersion matrix to use is dependent on 1) the format of the PCA problem at hand and 2) the goals of the analysis. If multiple fields contribute to the input data set (e.g., in a combined EOF) and have different units or levels of variability, some correction must be applied to the data to avoid having the field with greatest variance dominate the results. Normalizing each of the variables is one of the more common correction techniques, and is equivalent to using a correlation dispersion matrix. The considerations are more subtle when a single field is being analyzed.
PCA can be performed on observed fields, departure (anomaly) fields or normalized departure fields, corresponding to use of a cross-product, covariance or correlation dispersion matrix, respectively. There has been considerable discussion in the literature over the relative merits of each (Craddock 1965, Kutzbach 1967). While all three have been used in climatological applications, generally the covariance or correlation matrix is used (see Molteni et al. (1983) for a cross-products dispersion matrix analysis). Perhaps surprisingly, there is no simple relationship between the EOF/PCs of covariance and correlation matrices. Results from two analyses using the different dispersion matrices may differ substantially if individual variances among the variables are highly unequal. Use of the covariance matrix tends to isolate patterns in the regions of exceptionally high variability. Further, the relative differences between variables are preserved allowing one to ascertain quantitatively the degree of similarity among variables. If one is interested mainly in spatial covariability, particularly in areas of relatively little activity, then the correlation matrix should be used. I use the covariance matrix in these analyses.
Following preparation of the data, the eigenvectors, and eigenvalues are extracted from the dispersion matrix. The principal components are then computed by projecting the eigenvectors on the original data matrix (i.e., using equation (2), and the eigenvalues computed by solving equation (9).
An important concern in PCA is determining how many of the EOF/PCs represent signal and how many are unwanted "noise" (statistical significance). A second concern arises over the statistical uniqueness of each eigenvector (EOF), particularly when the number of samples (time observations) is small relative to the number of variables (grid points). Both issues are addressed through analysis of the eigenvalues. As noted earlier, the eigenvalues give the amount of variance explained by the transformed variables. The fraction of variance accounted for by a PC/EOF mode is calculated by dividing its eigenvalue by the sum ("trace") of the eigenvalues. If a covariance dispersion matrix was used, then the sum of the eigenvalues is equal to the variance of the dataset. If a correlation dispersion matrix was used, the sum equals the number of variables (grid points).
A large number of selection (or retention) rules to determine statistical significance of EOF/PC pairs have been derived, most based on the relative magnitudes of the eigenvalues. The number of EOF/PC modes to retain, and for that matter the meaning of the word "retain", is context-dependent. For example, if the intent of the PCA is to reconstruct a filtered field, retention of too few factors ("underfactoring") can result in loss of signal. The effects of "overfactoring" (i.e., retention of too many EOF/PCs) are increased computation time and retention of noise (underfiltering). If the goal is to physically interpret the EOF patterns, then one wants only to retain those that represent more than sampling error or environmental noise. In the matter of rotation, one interprets a select number of modes, but the number of EOF/PC modes retained for the rotation procedure is generally much larger. When rotation is applied to an underfactored solution, the discrete patterns blend together and interpretation of the results may not be possible (Richman 1981). Additional consideration of statistical significance comes from the eigenvalue separation (discussed below). Richman (1981) advises slightly overfactoring when the number of EOF/PCs to retain is indeterminate. Four of the more common selection rules are summarized below.
1) Scree test
One of the simplest tests, and probably the most widely used, is the scree test (Cattell 1966), which derives from a visual inspection of the eigenvalues plotted against their root number. A typical eigenvalue plot shows a steep slope over the first few roots and a gradual trailing off for the rest of the roots (the scree). The term scree derives from the resemblance to the rubble that forms at the foot of a mountain. Cattell hypothesized the scree represented unwanted noise and that only the EOF/PCs prior to the scree should be retained for further use. Identification of the exact root where the scree begins to form is rather subjective, however.
2) LEV graph
A technique similar to the scree plot is the LEV graph (Craddock and Flood 1969) where the Log EigenValue is plotted against its root. Experimental results indicate that a PCA of random data results in a straight line LEV graph (Farmer 1971). To determine the number of EOF/PC modes to retain, a straight line is drawn through the higher number roots, and those lying above the line are retained.
3) Guttman criterion
When the input dataset is normalized (or the correlation dispersion matrix used), the Guttman criterion (Guttman 1954) may be used. Under the Guttman criterion, only those EOF/PCs whose eigenvalue is greater than or equal to one is retained. The reasoning is that a typical normalized time series contributes one unit of total variance (l = 1.0), thus modes whose eigenvalue is greater than 1 represents a signal, the remainder noise. Horel (1981) found that rotation results were relatively insensitive to the number of modes retained as long as it was at least 25% of the number suggested by the Guttman criterion.
4) Rule N
A more recent test, based on a Monte Carlo procedure, is the Rule N test (Preisendorfer 1988). In this test, one creates a number of datasets (30-50 is typical) of random normal deviates with variances equal to those observed in the field under study. PCA is then conducted on these datasets and 95% confidence limits determined for each eigenvalue. The eigenvalues from the observed field analysis lying outside the 95% C.I. are then retained. It has been shown, however, that Rule N tends to underestimate the number of significant EOF/PC pairs (Preisendorfer 1988), possibly due to underestimation of the autocorrelation of the leading EOFs.
The uniqueness of EOFs is determined by the separation of their respective eigenvalues. North et al. (1982b) showed that when neighboring eigenvalues are similar in magnitude, an "effective degeneracy" occurs in which two or more EOFs become mixed. This degeneracy derives from the fact that any linear combination of the possible EOFs can also be an EOF. The result is that wildly different patterns will emerge from one sample to the next. Fortunately, there is both a diagnosis and solution to the problem. A useful rule of thumb was derived by North et al. (1982b) to determine whether neighboring eigenvalues might be subject to mixing. This rule says that an eigenvalue must be separated from an adjacent eigenvalue (Dab = la - lb) by more than its sampling error (dla=la(2/n)½); else they form part of a "degenerate multiplet" and their EOFs are a random mixture of the true EOFs.
The amount of separation between eigenvalues required to produce independent EOFs depends on the degrees of freedom, which in turn is inversely proportional to the number of cases n. Many climatological fields are autocorrelated in time, i.e., not statistically independent, thus the degrees of freedom is potentially far less than the number of observations. In such a situation, an "effective" number of degrees of freedom should be calculated and used in this rule of thumb. One commonly employed formula under the assumption of a red noise process (Laurmann and Gates 1977) is:
(11) 
where r is the lag 1 autocorrelation coefficient computed as:
(12) 
Most of the North Pacific spring SST field shows a lag 1 autocorrelation (r) of between .2 and .4 (Figure 1.3), thus resulting in a 33-57% decrease in degrees of freedom. In calculating the eigenvalue sampling errors for spring SST, I used a r value of 0.3, reducing the effective number of independent years from 43 to 23. Further discussion of this and other problems associated with autocorrelation is continued in the time series analysis section.
As the eigenvalue plot reaches the scree portion, none of the later eigenvalues is statistically separated from any of the others. This is another argument for not retaining higher numbered EOF/PC modes. It may happen that a series of degenerate multiplets occurs prior to, and continuing into, the scree. In such a situation, the truncation point for retention should take place between, not within, multiplets.
A useful method then of determining how many EOF/PC modes to retain is to combine a number of the tests listed above. I have found it statistically and visually informative to use the scree test, Rule N, and the North et al. (1982b) rule of thumb. These three tests can be plotted together, and the results for spring SST are illustrated in Figure 1.4. The Rule N confidence intervals were derived from PCA conducted on 100 random data sets of 23 years duration. The effect of autocorrelation in the datasets was accounted for by using a reduced number of years (23) than are in the actual dataset (43), based on the r value of 0.3 cited above. The random deviates were drawn from normal distributions with variances equal to that observed at the corresponding grid point in the spring SST data, and weighted by the square of the cosine of the latitude to offset the latitude effect.
As is often the case, the exact number of significant modes is somewhat indeterminate. The first mode, accounting for 39% (11.56/30.0) of the total field variance, is well separated from the other modes and much greater than the Monte Carlo 95% C.I. value. The eigenvalues for modes 2 and 3 are both greater than their corresponding Monte Carlo 95% C.I. value, but have significantly overlapping error bars. Mode 4 and 5 eigenvalues lie just below the 5%-95% C.I. bands and are almost equal in value. Modes 4 and 5 form a degenerate multiplet and must together be retained or rejected. There is a somewhat distinct break after the 5th mode with all higher modes having overlapping error bars and the scree having clearly begun. The main decision then is whether to interpret the first 3 after rotation. In this situation, I elected to retain only the first three modes. The degree of overlap between modes 2/3 and 4/5 is relatively minor and modes 4 and 5 are considered statistically insignificant by the Monte Carlo simulation. The first three modes of spring SST account for 63% of the total field variance.
If the input data are weighted (such as by the square root of the cosine) prior to formation of the covariance matrix, the resultant EOF values must be multiplied by the inverse of the weights to return the patterns to their actual relative values. In the traditional terminology and notation (but not that used here), the EOFs are nondimensional having been arbitrarily scaled to unit length. The PCs have the same units as the input data, i.e. raw or standardized values depending on whether the covariance or correlation dispersion matrix was used. Thus, each of the r EOFs is of unit magnitude, and information concerning the relative strength of the patterns must be derived from examination of the corresponding PCs, which in S-mode analysis are time series.
A potentially more useful display method, and that which I use in all analyses, is to weight the EOFs by the square root of their corresponding eigenvalue (the EOFs are then sometimes termed "PC loadings"). The resulting patterns are unchanged in shape from the original EOFs but show pattern strength, in actual units, corresponding to a PC value of +1 standard deviation. The first few EOFs will have a much larger range of values compared to later EOFs, reflecting their greater importance to the overall variability of the field. If the covariance matrix was used, the PC time series should be standardized, by dividing by the square root of its eigenvalue, to allow quick visual calculation of the magnitude of a perturbation pattern at any given time interval (i.e., EOF grid point value × PC standardized value = estimated deviation value due to that mode). With use of a correlation matrix, the PC is already standardized, and the resultant EOF × PC product is a unit deviate.
The strength of an EOF pattern is determined by the magnitude of both the PC and EOF. When displayed, the sign of the EOF is arbitrary. For instance, the first EOF of many fields tends to be a bullseye pattern all of the same sign (indicating a positive or negative departure of the field over the entire domain). When the EOF is displayed as positive, it means that the pattern is positive in value when the corresponding PC is also positive. A negative PC value means that the EOF is actually negative over the entire domain for that time interval. If one prefers, instead, to display an EOF as negative (say to match a pattern published in the literature), it is only necessary to reverse the sign of the PC as well.
The first three spring SST EOF patterns are illustrated in Figure 1.5 and the corresponding PC time series in Figure 1.6. The fraction of total field variance accounted for by each of the modes is shown with the EOFs. The interpretation of the spring SST EOFs is deferred to Chapter 2. It is common to also compute the correlation field between the PC time series and the observed time series at each grid point. Regions having a correlation greater than 0.5 are shaded in the EOF maps. These shaded regions indicate areas where the particular PC/EOF mode is highly indicative of actual variability at that point. This type of information can be extremely useful in selecting stations to use as indicators of a particular phenomenon, such as is done to extend a time series further back in time. Time series of SST extending before 1950 are available for several coastal Alaska cities (e.g., Juneau, Kodiak, Adak, King Salmon). On the basis of the PC correlations, the best SST measure of the physical process depicted in EOF1 would be from a northern Gulf of Alaska location such as Valdez or Cordova. Conversely, locations out along the Aleutian chain (Cold Bay, Dutch Harbor, Adak) would make poor proxies for this pattern of variability.
Because the PC is a quantitative measure of the temporal variability of the EOF, it is often subjected to further statistical tests. In particular, it is of interest to determine whether there is evidence of non-random behavior in the time series. For example, the PC can be used to search for cyclicity or regime-like behavior of a particular pattern, or maximum values of the PCs can sometimes be used to pinpoint periods of extreme departure from long term means. In this thesis, I apply time series analysis to many of the climatological PCs. An example of a time series analysis of a PC is given in the Time Series analysis section below.
All of the early, and most of the more recent, applications of the eigentechniques are limited to the initial EOF/PC solution. Beginning in the late 1970's, several researchers began to question the validity of such studies, particularly when the goal was physical interpretation of the eigenvector maps (Buell 1975, 1979; Horel 1981; Richman 1981). The controversy revolves around the issue of rotation.
Rotation involves taking a linear combination of the leading principal components (i.e., those retained according to the selection rules listed above). The purpose of rotation is to achieve what Thurstone (1947) termed "simple structure". The effect of a rotation to simple structure is to cluster within each mode a small number of high valued variables and a larger number of near-zero valued variables. The idea is that, within each mode, the high valued variables represent one "causal or summary influence" that does not affect the other (near-zero valued) variables. The orthogonality constraint of unrotated eigenvectors prevents a clean isolation of local (i.e., non-global) modes of variability.
There are two classes of rotation used in RPCA, each with many variants: orthogonal and oblique. In an orthogonal rotation, the PC coefficients remain orthogonal, however the eigenvectors are no longer necessarily uncorrelated. In an oblique rotation, the PCs are also allowed to be correlated. A number of climatological studies have found only minor differences in results when comparing techniques within a class (White et al. 1991, Richman 1986) or even between classes (Barnston and Livezey 1987). It is for this reason that the overwhelming majority of orthogonal analyses use the varimax scheme (Kaiser 1958) and oblique analyses (themselves only a fraction of the number of orthogonal analyses) use Harris-Kaiser (Harris and Kaiser 1964) or Procrustes (Hurley and Cattell 1962).
Debate over when and how to rotate EOFs and PCs continues to generate, at times, vociferous discussion (see, for example, the extended discussions of Mather (1971, 1972) and Davies (1971a, b, 1972); Richman (1986, 1987) and Jolliffe (1987); Richman (1986, 1993) and Legates (1991, 1993)). Perhaps the crucial point of these debates is that one should at least explore the possibility of examining alternative (i.e., rotated) solutions. I briefly summarize the arguments for (and against) rotation and then compare the results of a varimax rotation with the unrotated results shown previously
I begin with a few of the arguments against rotation and cite situations where rotation may be unnecessary. Mather (1972) claimed that "Rotations are irrelevant in principal component analysis", and should be reserved for factor analysis. This point of view derives from some researchers regard of PCA as "an empirical method for the reduction of a large body of data so that a maximum of the variance is extracted" (Harman 1967). Legates and Willmott (1983) concur with Mather (1972) stating: "Rotation should never be employed in PCA since its only effect is to destroy the maximum-variance explained property of the components." Legates (1983) later concluded that "If, however, a researcher is interested in extracting causes or modes of the variation ... then rotation is required because the goal is commensurate with factor analysis" (italics in original).
Another criticism of RPCA is the potential arbitrariness of the results due to the analyst's choice of 1) number of modes rotated and 2) type of rotation scheme used. Jolliffe (1987) termed the first problem "stability of rotated solutions", and it is related to the problem of sampling errors discussed earlier, and is revisited below. The results of a RPCA can change radically depending on the number of modes retained, whereas the PCA modes are invariant. It is true that a radical change in EOF patterns is noticeable when comparing the results of, say, a 1-mode retained versus a 2-mode retained analysis. However, this is likely due to the addition of significant variation due to an underfactored solution. Barnston and Livezey (1987) and Richman and Lamb (1985) found that the rotated solutions changed only marginally with retention of additional modes as long as all those representing signal had already been retained.
There are at least three situations in which rotation is unnecessary. One of the major uses of PCA is simply data reduction wherein one wishes to preserve the variability of a large set while retaining a much smaller number of variables. A second use of PCA, related to data reduction, is as a pre-processor for other multivariate techniques such as singular value decomposition, and canonical correlation, cluster and discriminant analyses. In this context, PCA serves to filter out unwanted noise and make the multivariate analyses less susceptible to sampling error from short time series (Preisendorfer 1988). Third, in situations where the first EOF/PC is dominant and the remaining modes of little explanatory power, physical interpretability of the EOF may be relatively straightforward (Richman 1993).
Summarizing a large body of evidence, Richman (1986) outlined four characteristics associated with unrotated solutions that "hamper their utility to isolate individual modes of variability": domain shape dependence, subdomain instability, sampling errors, and faithfulness to relationships embedded in the dispersion matrix.
Domain Shape Dependence
The problem of domain shape dependence was first identified by Buell (1975, 1979) and results from the orthogonality constraint of the eigenvectors. Buell (1975, 1979) demonstrated that eigenvectors from an analysis of a particular domain followed a predictable sequence (subsequently termed "Buell patterns") "which do not reflect the underlying covariation" (Richman 1986). For example, when the domain is rectangular/square and the first EOF is a "bullseye" pattern of one sign throughout the domain (+), the second EOF will be bipolar (+ -) with the zero line passing through the center of the EOF1 bullseye. Subsequent EOFs will be of shapes such as +/-, - + -, -/+, with more complicated shapes arising in higher order modes. Richman (1986) presents a table of 22 climatological PC analyses, of various meteorological parameters, that inadvertently illustrate the Buell pattern sequence. He concurs with Buell's (1979) contention that the unrotated results "should be looked upon with suspicion." Numerous other authors, noting the appearance of Buell-like patterns in their unrotated results have similarly recommended rotation before interpretation (e.g., Cohen 1983, Molteni et al. 1983, Lins 1985, Crane and Barry 1988, Nicholson and Nyezi 1990)
Close examination of the first three spring SST EOFs (Figure 1.5) shows some evidence of the classic Buell patterns. The first EOF is a dipole pattern with a zero line following the coast. The second EOF is a mainly positive (bullseye) pattern with its maximum located over the EOF 1 zero line. The third EOF is also bipolar but with its zero line situated perpendicular to the zero line of EOF 1. Two factors likely cloud the Buell patterns in this example. First, the domain used here is not square or rectangular due to the absence of Bering and Okhotsk Sea data. Secondly, the eigenvalues for the second to fifth EOFs overlap resulting in sampling degeneracy (North et al. 1982), and mixing of the true EOF patterns.
The appearance of Buell patterns does not automatically mean that the results are just mathematical artifacts. One can easily visualize recurring natural phenomena in which the true variability matches a Buell pattern (e.g., opposing high and low pressure cells). In such a situation, rotation would preserve the patterns if they displayed simple structure in a plot of their unrotated variables. A rotated solution has been posited as being free of the domain shape influence (Kaiser 1958, Horel 1981, 1984).
Subdomain Instability
A corollary to domain shape dependence is the problem of subdomain instability (Horel 1981, Richman 1986). Ideally, the dominant EOF patterns, identifying coherent regions of variability, should retain their shape as the study region is varied in size. For example, a pattern identified for the North Pacific should also be evident if the analysis is then conducted on the Northern Hemisphere. In an intriguing study, Richman and Lamb (1985, 1987) conducted a PCA on rainfall data for the central United States. Buell patterns were found for the first ten unrotated EOFs, thus the first EOF showed a pattern of positive rainfall across the entire central US; the second EOF showed a contrast between the northern and southern halves, the third between the east and west halves, etc. Such patterns might be real, and if so, should be evident in an analysis of subdomains of the region. They next split the region into northern and southern halves and conducted separate analyses on each half. What they found was that the same sequence of Buell patterns was "simply transferred and compressed into the subareas". Centers of activity identified in the full domain PCA were not found in the separate domain analyses. In a varimax rotated PCA of the same three domains, the problem of subdomain instability was eliminated as the separate domain patterns matched those identified in the full domain analysis.
Sampling Errors
A number of authors have identified problems with unrotated EOFs when neighboring eigenvalues are closely spaced (North et al. 1982b, Storch and Hannoschöck 1985). If a multiplet (EOFs whose eigenvalue error bars overlap) is present in the retained EOF/PCs, then it is likely that the resultant patterns represent a mixture of the true patterns and no physical interpretation should be made of any of the EOF patterns (Storch and Hannoschöck 1985). Richman (1993) noted that in some studies in which the classic Buell patterns didn't arise, it may have been due to mixing of the patterns that resulted in new ones that were less recognizable. Rotation of EOFs has been found to result in stable, interpretable patterns (e.g., Hsu and Wallace 1985) due to the fact that the rotated EOFs have much lower sampling errors than their unrotated counterparts (Horel 1984, Richman 1986, Cheng et al. 1995).
Faithfulness to Relationships Embedded in the Dispersion Matrix
The purpose of PCA is to extract patterns of variability that reflect known natural occurring phenomena. Another aspect of unrotated solutions is that they are dominated by the maximum variance explained requirement. This requirement applies to the entire domain and is one of the reasons the bullseye pattern arises as EOF1 in many climatological studies. In a situation where several local areas within a large region have patterns of variability, the result of PCA will be to merge these unique signals into a single signal with no physical interpretability (Karl et al. 1982). By rotating, the global explanation requirement is eased and local features are more easily identified. A convenient means of testing the faithfulness of the first few EOFs is to compare them with anomaly maps of the time periods identified by the PCs as being strongest. One can also calculate correlation maps between EOFs and the original data to determine if individual areas of variation are being identified.
The idea behind a varimax rotation of principal components is to simplify the eigenvectors by concentrating the variability within each eigenvector to a small fraction of the variables. Ideally, each eigenvector would be composed of a few high valued variables, while the remaining variables would be near zero-valued. Further, each variable would have a high loading on one eigenvector and low loadings on the others. The criterion for accomplishing this is the statistic
By maximizing the variance of a fourth power sum, the coefficients tend towards 0 or 1. Maximizing this criterion is an iterative procedure that works on all possible pairs of retained components. Of all rotation criteria, varimax is the most widely implemented. In this study I use a FORTRAN version originally written by Kaiser (1958 and the implementation in the SPSS (1993) Professional Statistics microcomputer package.
On the basis of eigenvalue analysis, I decided to display (as statistically significant) the first three EOF/PC modes of spring SST after rotation. However, the first 10 PCs were rotated to prevent stability problems associated with "underfactoring" (Richman 1981). The rotated empirical orthogonal functions (REOFs) are illustrated in the right side of Figure 1.5. The most striking change is the splitting of EOF 1 into two components REOF 1 and REOF 2. The effect of rotation was to isolate two apparently separate modes of variability that were "smeared" together in unrotated PCA. This is another illustration of the tendency towards unrotated global (i.e., having a pattern with structure over the entire analysis domain) solutions that occurs in most climatological PCA applications and why one should explore alternative solutions. REOF 3 also bears some resemblance to EOF 2. The other effect of rotation is to redistribute the amount of variance among the modes retained in the analysis. While the cumulative variance for all 10 modes is the same for rotated and unrotated solutions, the effect of rotation is to distribute the variance more evenly among the modes. The corresponding rotated principal components (RPCs) are illustrated in the right side of Figure 1.6. It is interesting to note that RPC 1 and RPC 2 each show a similar time history to unrotated PC 1: predominantly negative until the mid 1970's and positive afterwards. The meaning of these results will be discussed in Chapter 2.
A principal feature of historical/retrospective analysis is that one of the measurement axes is time. Intrinsically, observations on a phenomenon (or process) over time are correlated with, or dependent on, their past. This has two important ramifications; first, the order of observations is important and second, the assumption that successive observations constitute independent samples is invalid. For example, although the t-test for equality of means is relatively robust to violation of the normality assumptions (Wei 1990), it has been shown to be invalid when there is even "moderate dependence between observations" (Box and Tiao 1965). A second example involves correlation between two time series. Katz (1988) has shown that internal autocorrelation leads to an apparent lead/lag relationship between two time series when none actually exists.
A set of statistical tools and techniques, collectively referred to as time series analysis, (TSA) has been developed to analyze data collected sequentially over time. There are three main goals of time series analysis. The first goal is to diagnose the nature of time-related behavior within and between time series. The detection of autocorrelation within a series (univariate analysis) is accomplished via the sample autocorrelation functions (ACF) and sample partial autocorrelation function (PACF). Time-lagged correlation between time series (multivariate analysis) is identified via the sample cross-correlation function (CCF). The second goal of TSA is to model the dynamics of a time series as a function either of its own past history (univariate TSA) and/or the history of other explanatory time series (multivariate TSA). The third goal of TSA is prediction of the future behavior of a time series given information concerning its present state. In this review, I restrict my coverage to the first two goals as I make no attempt to formally predict the future of either climate or Alaskan salmonid populations.
In this dissertation, I fit AutoRegressive Integrated Moving Average (ARIMA) and Intervention models to a wide variety of physical and biological time series. These models are first used to explore the nature of internally generated variability (AR and MA processes) and response to external forcing (interventions). The models then serve as prewhiteners in the ensuing cross-correlation analyses. In this chapter, I provide a brief outline of several important techniques (and associated notation) from the large field of time series analysis which I then demonstrate by way of example. The subjects I cover are correlation (auto-, partial auto-, and cross-), ARIMA and Intervention modeling, and prewhitening. Those seeking a more theoretical development should consult one of the numerous texts available including the seminal works on ARIMA model formulation (Box et al. 1994) and intervention analysis (Box and Tiao 1975).
In most cases, we have only a single realization (time series) of a process. Without having the ensemble of possible realizations of the process generating the time series we must rely upon sample estimators to characterize the properties of our time series. A time series is characterized by its mean m, its variance g(t), and its autocorrelation structure. In order to develop inferential procedures in TSA, a time series must be stationary: i.e., the mean and variance must be constant over time.
A time series is further said to be weakly, second order covariance stationary if:
Given these conditions, we can proceed with our single realization since the relationship between points in the series depends only on their relative, and not absolute, position in time. In other words, a different realization of the process, collected during some other time interval would have the same correlation structure. We use the following estimators for the mean and autocovariance function (which is the variance at lag k = 0):
If a time series is nonstationary in either the mean or variance, simple procedures exist to transform the time series to stationarity. To stabilize the mean of a series, it is common to difference the data:
where zt is the series resulting from differencing yt at lag k and results in a series of length n-k. The most common forms of differencing are lag 1 ("first") differences which are applied to "random walk" time series, i.e., series with a stochastic trend; and "seasonal" differences for data containing a regular cycle (e.g., hour or season). A series can be differenced multiple times, such as to remove an upward trend and an annual cycle. For series containing a dth-degree polynomial trend, applying first differences d times will remove it.
Variance non-stationarity occurs when the variance of a time series changes with time or with level (e.g. greater variance at high values). Box and Cox (1964) introduced a class of variance stabilizing transformations:
Typically, l is restricted to the range -1 £ l £ +1. In practice, for fisheries purposes, l values of 0, .5, and 1 are all that are considered. A value of zero, corresponds to a natural logarithm transformation and is appropriate when the variance increases quadratically with the mean. If the variance appears to increase linearly with the mean, then the square root transformation (l = 0.5) is appropriate. A l value of 1 is equivalent to no transformation. While Box and Cox (1964) use a maximum likelihood estimation procedure to estimate l, a visual examination of the potential transformations is usually adequate (Kendall et al. 1983).
When transformation of the original time series is deemed necessary, the variance stabilization is first applied followed by any requisite differencing. The Egegik sockeye salmon run time series analyzed in Chapter 3 was found to require both types of stabilizations. In Figure 1.7, I plot the original time series, the ln transformation and the differenced (of the ln transformation) time series. It is evident that both transformations at least partially accomplished their tasks. The ln transform reduced the much larger variance in the second half of the time series. The ln-transformed series is characterized by a negative trend in the first half of the series and a positive trend in the later half. First differencing removed both trends. An important point to make here is that a time series containing two distinct levels is not considered to be non-stationary in the traditional sense. An abrupt change in level is detected and modeled using Intervention Analysis, discussed below.
Once a series has been properly transformed, the next step is to characterize its autocorrelation structure. Two tools are used for this purpose, the autocorrelation function (ACF) and partial autocorrelation function (PACF). In a TSA, we are working with a finite length realization of a process, hence we can only obtain an estimate of the autocorrelation. As a result, we use the sample ACF and sample PACF. The most commonly used version of the sample ACF is:
It can be seen the ACF is restricted to the range -1 £ 
£ 1, and that the function is symmetrical
about lag 0. It will be of interest to determine whether the values of
the ACF are significantly different than zero at lags other than 0. To
do so requires an estimate of the variance of the ACF, the sample approximation
of which was shown by Bartlett (1946) to be:
A plot of the autocorrelations against lag (the "ACF plot", sometimes termed the sample correlogram) generally includes plots of ± 2 standard errors, i.e. two times the square root of Eq. (21). Due to its symmetry, and its perfect correlation at lag 0, An ACF plot usually includes only positive lags and is limited to £ n/4.
When there is high autocorrelation within a series, there will be a "carryover effect" at neighboring lags making it appear that there is significant correlation at several lags. While such a relationship is possible, it may also be spurious. The sample partial autocorrelation function is used to estimate the correlation between zt and zt-k after removing the mutual linear dependency on the intervening values zt-1, ..., zt-k+1. There is no single formula for the sample PACF, as it is lag-dependent. By definition, however, it is the following ratio:
where 
is the estimate of the last coefficient in an AR(k) model.
The standard errors of the PACF estimates are generally calculated using the following approximation, which is based on an underlying assumption of a white noise process:
The interpretation of ACF and PACF plots can be as much art as science. As is illustrated in Chapter 3, natural processes can give rise to unusual plots, indicative of complex dynamics. Still, the first step in a TSA is examination of the ACF and PACF plots as often the nature of the underlying process can be deduced from certain patterns. For example, a time series that is uncorrelated in time is referred to as a white noise time series and is characterized by ACF values all contained within the ± 2 standard error bars. An ACF that alternates between positive and negative values indicates an anti-correlated process. A non stationary time series is characterized by a ACF in which the values die out very slowly. ACF and PACF plots for theoretical ARIMA processes are provided in most standard time series texts, though most warn that plots from real data are never as clear.
One of the most common ACF/PACF forms seen in environmental modeling is one in which both the ACF and PACF are significant only at lag 1. This possibly indicates a "red noise" process, i.e., one that slowly varies in time such as a low frequency cycle. An alternate explanation for such a combination of ACF/PACF plots is given in the section on prewhitening. The ACF and PACF plots for the spring SST 1st principal component (Figure 1.8) show such a pattern. The lag 1 coefficient (which is the same for the ACF and PACF) is .480, which is highly significant compared to the standard error value of .152. The lag 2 and lag 3 coefficients both appear to be nearly significant as they are approximately twice the value of their respective standard errors. The PACF, however, shows that the high correlation at lags 2 and 3 results from the "carryover" effect noted above. The traditional interpretation of these ACF/PACF plots is that the dominant mode of spring SST variability is significantly driven by its immediate past history. One way of looking at this is to say that SSTs are anomalously high this year because they were last year and, therefore, they will be high next year.
The cross-correlation function (CCF) is used to detect the presence of a lead/lag relationship between two time series, denoted xt and yt. The same stationarity requirements apply to each series as applied to univariate series prior to calculation of the ACF and PACF. Similarly, it is useful to regard the two time series as realizations of a bivariate stochastic process. We say they are jointly stationary if the cross-covariance between xt and yt is a function only of the time lag k.
The cross-correlation function is defined as
and the sample cross-correlation function is estimated by:
where
To test whether there exists significant cross-correlation, we compare the sample CCF lag values with their standard errors. Bartlett (1955) derived a formula to compute the covariance between two cross correlation estimates and, hence, the variance and standard errors at individual lags. The formula is complicated but, under the hypothesis that xt and yt are uncorrelated and the xt series is white noise, the variance can be approximated by:
An important aspect of the CCF is that, unlike the ACF
and PACF, it is not symmetric about the origin, i.e., 
.
A CCF plot (often called a cross-correlogram) includes both positive and
negative lags, thus indicating both the strength and direction of an association.
When the CCF is computed as 
,
i.e., correlating the series xt with lagged values of yt,
positive lags on the CCF plot show x leading y, negative lags show y leading
x.
In practice, xt is not white noise and the resultant CCF can be highly distorted, usually indicating significant lagged correlation which is actually spurious. This is a highly significant point often ignored in fisheries oceanography work on climatic influences on biological processes. A more detailed examination of this problem and the necessary steps to remove the spurious correlation ( a process known as "prewhitening") is discussed in detail in Section 2.3.3. Examples of CCFs are also postponed until the discussion on prewhitening.
The use of time series analysis to model fish population dynamics has increased in recent years. Most of the theoretical development and initial application has taken place in the econometric and business forecasting literature. Recognition of the potential applicability to ecological problems appears to have begun with Moran (1949).
There are five classes of commonly applied time series models (Jenkins 1979). The simplest, and most widely known, comprise the so-called Box-Jenkins ARIMA univariate models. Simple ARIMA models utilize only the history of the time series to "explain" its observed variability. The second class comprises the transfer-function noise (TFN) models, which relate an output-series variability to both its own history and that of one or more explanatory variables. A third class, related to TFN models, comprises intervention models which incorporate the effects of unusual events, natural or human-made, to modify ARIMA models. The other two classes comprise multivariate models. Multivariate stochastic models permit feedback among several time series and are often referred to as vector ARIMA models. The final class includes explanatory variables giving a multiple input/multiple output mode and are sometimes referred to as multivariate transfer-function models.
In addition to these time series models, there has been a parallel development of frequency-domain models, principally in the engineering literature. In the frequency-domain models, processes are modeled as combinations of cosine waves. While theoretically translatable to time-domain models, there have been few applications in ecology. More recently, state-space models have generated a great deal of attention. In state-space, or more generally, structural modeling, a time series is decomposed into linear, seasonal, and irregular components (Harvey 1989). The central feature of structural models is the use of the Kalman filter (Kalman 1960; Kalman and Bucy 1961) for parameter estimation and forecasting. The principal difference between traditional time series and structural models is the manner in which the error component is modeled. Though neither method has emerged as clearly superior, structural models are likely to receive increased attention.
The first published use of time series modeling in the fisheries literature was Dunn and Murphy (1976) and Murphy and Dunn (1977), who used univariate and transfer-function models to forecast fish catch in an Arkansas reservoir. Univariate and/or transfer-function models have been used to model the population dynamics of American lobster (Homarus americanus; Boudreault et al. 1977, Fogarty 1988a, Campbell et al. 1991), rock lobster (Jasus edwardsii; Saila et al. 1980), skipjack tuna (Katsuwonus pelamis; Mendelssohn 1981), yellowtail flounder (Limanda ferruginea: Kirkley et al. 1982), menhaden (Brevoortia patronus; Jensen 1985), haddock (Melanogrammus aeglefinus; Pennington 1985), Alaskan salmon (Quinn and Marshall 1989; Noakes et al. 1987), winter flounder (Pseudopleuronectes americanus; Jeffries et al. 1989), blue whiting (Micromesistius poutassou; Calderon-Aguilera 1991), pilchard (Sardina pilchardus; Stergiou 1989), and striped bass (Morone saxatilis; Tsai and Chai 1992). Intervention analysis has been applied to Dungeness crab (Cancer magister; Noakes 1986), geoduck clams (Panope abrupta; Noakes and Campbell 1992), power plant impact on yellow perch (Perca flavescens) and alewife (Alosa pseudoharengus; Madenjian et al. 1986), and to forecast invertebrate yield (Fogarty 1988b). Vector ARIMA models have been applied to Great Lakes pelagic species (Cohen and Stone 1987; Stone and Cohen 1990) and multivariate transfer function models were used by Mendelssohn and Cury (1987, 1989) to explore catch per unit of effort in Ivory Coast pelagic fisheries.
Use of formal time series methods in climatological research has also been somewhat limited. The theory of stochastic climate models (Hasselmann 1976, Frankignoul and Hasselmann 1977, Lemke 1977) is based on a random walk model with feedback where climate acts as a long period integrator of short-term random weather excitation. As noted earlier, recent developments in modeling of dynamic climate processes have merged multivariate PCA and ARIMA techniques to form PIP, POP, and PPP models (Hasselmann 1988, Bürger 1993). Time series modeling has been used to analyze rainfall (Chang et al. 1984), global warming (Woodward and Gray 1993), nonstationary climate processes (Rao et al. 1992), the greenhouse effect (Tol and de Vos 1993) and to forecast wind speed (Brown et al. 1984). Statistical problems that arise in climatological research, and the need to employ a time series approach, has been addressed by Katz (1988a, 1988b, 1992), Brown and Katz (1991), and Zwiers (1990).
ARIMA and intervention models have several different representations. We employ the following notation:
(31) 
is the discrete
time series, which may be transformed to stabilize the variance using the
Box-Cox (1964) power transformation. The most common transformations are
square root (l=0.5), natural logarithm (l=0.0),
and inverse (l=-1.0). No transformation is equivalent
to a lambda value of 1.0. If required, a power transformation must be done
as the first step in time-series modeling.
is an "integrating
factor" (the "I" in ARIMA), better defined as a differencing
operation to induce stationarity in the mean of a series. The number of
differences taken (which can be at various lags) is indicated by d.
If required, differencing is the second step in ARIMA modeling.
is a seasonal
integrating factor(s) where s is the lag at which the Dth
seasonal difference is taken. While seasonal models are generally applied
to weekly, monthly, quarterly, etc. data, they may also be applied to non-seasonal
data that exhibit seasonal (i.e., periodic) behavior.
plays different
roles depending on the value of d (order of differencing). For d
= 0, q0 is equal to the estimated
mean of the transformed input series multiplied by the sum of the autoregressive
components and moved to the right-hand side of the equality. For d
³ 1, q0
is called the deterministic trend and is often omitted unless clearly called
for (Wei 1990, p. 72).
at is a random error component assumed to be normally independently distributed with mean 0 and constant variance s2a.
B is the backshift operator. By convention it is a special notation used to simplify the representation of lagged values: Byt = yt-1, Bsyt = yt-s. Note also the following definition: Ñ = 1 - B, thus differencing is often represented by: Ñyt = (1 - B)yt.
is the autoregressive
polynomial of the form (1 - f1B -
f2B2 - ... - fpBp).
The term "autoregressive" is in reference to how the value of
y is being regressed on its own past values plus a random shock,
thus relating the present value of a process to a linear combination of
its past values. An autoregressive process can be written as yt
= f1yt-1 + f2yt-2
+ ... + fPyt-P
+ at.. An autoregressive process of order p is abbreviated
AR(p), and lower orders than p need not be non-zero.
is the multiplicative
seasonal autoregressive polynomial of the same form as the non-seasonal
polynomial. Multiple seasonal autoregressive components may be included
in the model, each of seasonality S. The subscript P identifies
the presence of a seasonal component, and all coefficients other than that
of the seasonal lag are set equal to 0.
is the moving
average polynomial of the form (1 - q1B
- q2B2 - ... - qqBq).
The moving average term models the persistence of random effects over time
and can be written as yt = at + q1at-1
+ q2at-2
+ ... + qpat-p.
A moving average process of order q is abbreviated MA(q), and lower orders
than q need not be non-zero.
is the multiplicative
seasonal moving average polynomial of the same form as the non-seasonal
polynomial. Multiple seasonal moving average components may be included
in the model, each of seasonality S. The subscript Q identifies
the presence of a seasonal component, and all coefficients other than that
of the seasonal lag are set equal to 0.
represents
the jth intervention and is analogous to a dummy variable in regression.
Interventions can be either step (I = 1 for t ³
T, I = 0 otherwise) or pulse (I = 1 for t = T, I = 0 otherwise)
functions. A step intervention indicates a permanent shift in the mean
of a series, while a pulse indicates a one-time shock. There are several
different system responses to step and impulse interventions, such as an
abrupt permanent step, a step decay, and impulse decay.
is a polynomial
of the form (w0 - w1B
- w2B2 - ... -
wsBs) representing
the initial impact of the intervention.
is a polynomial
of the form (1 - d1B - d2B2
- ... - drBr) representing
the long-term impact of the intervention.
models the
delay in response associated with a particular intervention.
Nonseasonal ARIMA models use the notation (p, d, q) to compactly represent autoregressive, difference, and moving average orders. Seasonal models are expressed as (p, d, q) x (P, D, Q)S, with each seasonal component separately represented. Thus, a (1, 0, 5) model indicates the presence of additive lag 1 AR and lag 5 MA terms with smaller lag MA terms possibly present. A (1, 0, 0) x (0, 0, 1)5 model also has lag 1 AR and lag 5 MA terms, but the parameters are multiplicative rather than additive.
Univariate time-series model building, in the methodology of Box and Jenkins (1976), proceeds in the following fashion:
1) Model Identification. In this step, tentative models are identified. Determination of the need for power transformation (for variance stabilization) and differencing (to render the series stationary in the mean) are first evaluated. Plots of the autocorrelation and partial autocorrelation functions (ACF and PACF respectively) of the possibly transformed series are examined to assist in determining the order of the AR and MA components (Box and Jenkins 1976). Several other identification tools are also available, such as the extended sample autocorrelation function (ESACF; Tsay and Tiao 1984), generalized partial autocorrelation coefficient (GPAC; Woodward and Gray 1981) and the prediction variance horizon (PVH; Parzen 1981).
2) Parameter estimation. Following selection of a potential model(s), estimates of the parameters are calculated. Access to time-series software is almost essential as ARIMA model parameters must be fitted using a nonlinear estimation routine (though the models themselves are usually linear). Maximum likelihood procedures, usually based on the Cholesky decomposition or the Kalman filter, have been developed as an alternative to the early methods of least squares and approximate likelihood utilized by Box and Jenkins (1976). Standard errors are also computed, and parameters judged to not be significantly different from zero can be dropped. The remaining parameters are then re-estimated.
3) Model diagnostic checking. With a tentative model selected and parameters estimated, the adequacy of the model must be assessed to determine if model assumptions are met. One basic assumption is that the residuals at form a white-noise series. A common test is the portmanteau test of Box and Pierce (1970), which uses the residual ACF to test the joint null hypothesis that all serial correlations are equal to zero. It is also common in time-series analysis that several models may be adequate in the sense that the model residuals are reduced to white noise. Several model selection criteria have been developed to assist in model selection. In this analysis, I compare competing models using five criteria: mean absolute error (MAE), which measures the average one-step-ahead prediction error; the unbiased residual variance s2a, equal to the residual sum of squares divided by degrees of freedom; the coefficient of determination r², which is the amount of variance "explained" by the model; Akaike's Information Criterion (AIC; Akaike 1974); and Schwarz's Bayesian Criterion (SBC; Schwarz 1978). The AIC and SBC are performance statistics that balance statistical fit with model parsimony. The SBC utilizes a larger penalty function than the AIC, thus often suggesting a model with fewer parameters. Formulas for the model diagnostic and selection criteria are contained in the Appendix.
As an example of the Box-Jenkins iterative fitting procedure I will use the spring SST PC1 time series. The ACF and PACF plots shown earlier each show that only the lag 1 coefficient is significant. On this basis, three ARIMA models are potential candidates: (1,0,0); (0,0,1) and (1,0,1). In addition to the AR and MA parameters, one must also determine whether a constant (intercept) is required. Therefore, six univariate models - designated U1 through U6 - were fit (Table 1.2). Models U1, U3, and U5 each contained a constant and in all three cases, the constant was not significantly different than zero. This is unsurprising as the time series is a standardized principal component (i.e., of mean zero and unit variance). Thus, the number of candidate model is reduced to three.
To choose among the three remaining models, the next step is to examine the performance statistics (Table 1.3). The residuals from all three model fits pass the portmanteau white noise test - a necessity for an acceptable model. Model U4 (the (0,0,1) model) is clearly outperformed by the two other models, having a higher mean absolute error, residual sum of squares, unbiased residual variance and much lower r2 value. Choosing among models U2 and U6 requires more careful examination. Model U6 (1,0,1) does provide a better fit than model U2 (1,0,0), however this result is to be expected given the extra parameter. To decide if the fit is sufficiently better to warrant retaining the more parameterized model, we examine the two fitting criterion based on the log-likelihood value. Here, U2 has both a lower AIC value and SBC value. I conclude, therefore, that the best model for this time series is the (1,0,0) with no constant. This is a simple lag 1 autoregressive model suggesting that spring SST does have a significant relationship with its immediate past history. I note, however, that this model "explains" only 20.6% of the total variance in the series as indicated by the r2 value.
The purpose of intervention analysis is determine the effect of extraordinary events on a time series of interest. Commonly cited examples of intervention analysis include studies on the effect of labor strikes, advertising promotions, power outages, and introduction of environmental regulations. These interventions represent a different dynamic within a time series that overlays the internal ARIMA structure. Indeed, it is the influence of the autoregressive and moving average components that make identification of the intervention effect difficult.
Intervention events may be either permanent or temporary and their influence can be gradual or abrupt. It is also possible for there to be both a temporary and permanent effect (such as large initial effect with more modest permanent effect). In all cases, the interventions occur at a single time step with variable response in the impacted time series. The permanent events are referred to as step interventions and result in a permanent change in the mean level of the process under study. Temporary events, termed pulse interventions, have a brief impact on the level of the process but do not cause a long-term change in the mean level. Examples of the variety of responses to pulse, step and pulse-step combined interventions are illustrated in Figure 1.9.
In intervention analysis, the correlation structure is initially assumed to be unaffected by the interventions that are modeled as deterministic functions of time. Once the best ARIMA model has been selected, the three-step modeling sequence is repeated to identify and test the significance of interventions. The original intervention methodology developed by Box and Tiao (1975) permitted estimation of intervention effects when the timing and form of the interventions was known a priori. The intervention parameters are estimated via maximum likelihood simultaneously with the ARIMA parameters and standard errors calculated to determine significance. This is the technique implemented in most TSA software packages.
To handle the situation where the number and timing of potential interventions are unknown, Chang and Tiao (1983) proposed an iterative detection technique using a likelihood ratio test. Interventions are identified in a stepwise fashion beginning with the residuals from the univariate model. Following detection and estimation of an intervention, model parameters are estimated and the resultant intervention model compared with the univariate model using the criteria cited above. The new model residuals can then be re-analyzed for evidence of other interventions.
It should be noted that testing for different types and timings of interventions increases the probability of identifying a spurious intervention. However, I minimize this risk in two ways. First, my use of the AIC and SBC performance statistics helps insure that the increased model fit resulting from an intervention is significant. Secondly, in all cases, I specify the timing of an intervention and then test for its significance rather than blindly searching for the time period when the model performance statistics are maximized. Two software packages, AUTOBOX (Automatic Forecasting Systems, Inc. 1992), and SPSS Trends (SPSS, Inc. 1993), were used for all analyses. A good general review of intervention models is contained in Wei (1990), while Noakes (1986) discusses the applicability of intervention analysis to fisheries problems.
To demonstrate intervention fitting, I continue the example of the spring SST PC1 time series. The best univariate model for this series was found to be an ARIMA(1,0,0) model. On the basis of the climate analysis of Chapter 2, I hypothesize that a step intervention should be incorporated into a model for this series and it should take effect in 1977. The first intervention model (IV1) is a (1,0,0) model with no constant and a step intervention in 1977. The intervention coefficient is significant for this model and it clearly outperforms the non-intervention models (Table 1.3). Because no constant is specified the intervention effect may be constrained since it can only change from 0 to the estimated mean following the 1977 intervention. Thus, in model IV2, I add an intercept parameter and refit the model. The addition of the intercept, and resultant increase in the estimated intervention effect, has the added effect of reducing the significance of the lag 1 autoregressive parameter. So the final model I fit, model IV 2, contains no autoregressive or moving average parameters (i.e., a (0,0,0) model), only an intercept and intervention term. While several of the statistics slightly favor model IV2, the principal criteria - the AIC and SBC - both strongly favor model IV3.
The interpretation of this result is rather interesting. What it says is that large scale variability in spring sea surface temperature is not best explained as a slowly varying autoregressive process. Instead, it suggests that 1) SST alternates between cold and warm regimes, 2) the transition from one regime to another is abrupt, and 3) variability is random (in the time domain) within a regime. Since 1977, SST has randomly varied around some mean level rather than slowly increasing and decreasing and this mean level is considerably warmer than the pre-1977 level. This topic is covered in much greater detail in Chapter 2.
The CCF is extensively used to test for the presence of a lead/lag relationship between time series. In particular, I use the CCF for two different types of tests. In the first situation, two time series are tested for a lead/lag or feedback relationship. Thus, either series is presumed capable of affecting the other, such as may be the case between environmental variables In other words, a situation in which it is not known which direction the influence acts. The four possible relationships between two time series, X and Y, and the form of the CCF are:
(i) X leads Y
and
(ii) Y leads X
and
(iii) Feedback exists between X and Y
and
(iv) Only a contemporaneous relationship exists between X and Y
and
The second situation is concerned with an explanatory (and, therefore, possibly predictive) relationship particularly between environmental and biological variables. In such a situation only cases i) and iv) above are possible or sensible.
As shall be illustrated below, the presence of autocorrelation within the individual time series can seriously contaminate the CCF, giving the appearance of a significant relationship when none actually exists. This phenomenon, though well documented in the time series and econometric literature (Sims 1972, Pierce and Haugh 1977), has only recently begun to be advertised in the environmental sciences (e.g., Katz 1988b, Brown and Katz 1991), and has received little attention in the fisheries and fisheries oceanography literature.
The principal means of dealing with this problem is via a technique termed "prewhitening" which, in essence, filters the time series for autocorrelation. The CCF is then computed from the residuals of the prewhitened series. It has not been uncommon to find in the literature (e.g., Niebauer 1984, Gordon 1986, Salmon 1992) instances of smoothing time series prior to computation of the CCF, the justification being that 'noise' might otherwise obscure a relationship. Smoothing, in fact, introduces autocorrelation into a time series and, thus, is just the reverse operation of what should be done.
The goals of this section are threefold. First, I briefly illustrate mathematically how autocorrelation within two time series can distort the CCF between those time series. Secondly, I demonstrate by example how the CCF between some sample physical/biological time series changes (or doesn't) with prewhitening, depending on the autocorrelation structure. Finally, I extend prewhitening to include removal of intervention effects using simulation to illustrate how uncorrelated series that undergo separate simultaneous interventions (such as occurs during regime shifts) will also bear the illusion of a significant relationship. One crucial point that should emerge from this discussion is that a true lead/lag - and, therefore, possible cause/effect - relationship between two time series will survive, indeed, be clarified, by the prewhitening process.
Two different forms of prewhitening exist, depending on which of the two situations described above (i.e., bidirectional or unidirectional relationship) one encounters. The justification for prewhitening is also slightly different in the two cases. They are separately described below.
Situation 1. Potential bidirectional relationship
To illustrate the effect of autocorrelation on the CCF in this situation, we define the following two AR(1) time series:
To simplify matters, assume that both X and Y have been
standardized, i.e., scaled to zero mean and unit variance. The two residual
series, 
and 
are white noise processes with mean zero and variances 
and 
respectively (see
Box et al. 1994), Chapter 2 for derivation), where 
and 
are the lag 1 autocorrelation coefficients. If the two white noise time
series are significantly related contemporaneously (i.e., lag 0 as in relationship
(iv) above), but not at any other lag, the CCF between the AR(1) series
X and Y can be shown to be (Bartlett 1955):
and
These formulae illustrate that the CCF is directly influenced
by the level of autocorrelation within the time series. In the presence
of substantial autocorrelation, the true relationship between X and Y,
whether lead,. lag, feedback, or contemporaneous, would be masked, particularly
when combined with sampling error. Furthermore, the contemporaneous lag
0 relationship between 
and 
is also distorted by the autocorrelation coefficients (Katz 1988b):
where 
is
the lag 0 cross correlation between the two noise processes.
To remove the internal autocorrelation within both time series, both series are fit to (filtered by) Box-Jenkins ARIMA models using the methods described earlier. The residuals from the time series models must be white noise and these are the series used in calculating the CCF. Importantly, in the bidirectional case, separate models are fit for the two series. This method of prewhitening is often termed "double prewhitening" (Wei 1990, p. 299) due to the use of separate filters. A thorough discussion of double prewhitening with more detailed examples can be found in Haugh (1976).
Situation 2. Unidirectional relationship.
The CCF is closely related to the impulse response function (IRF) used in construction of multivariate transfer function models, relating a response variable to past and present values of an explanatory variable. A transfer function model can be written as follows (Wei 1990):
where 
is the IRF and n is a noise process. The CCF and IRF are related
as follows:
Here, it can be seen that the autocorrelation structure
of the forcing variable X contaminates the relationship. If X were a white
noise series (i.e., 
=
0 for all k), the CCF would then be directly proportional to the IRF.
In the unidirectional situation, a time series model/filter is developed for the explanatory variable and both series are filtered using this model. This is known as simple prewhitening, and is the more traditional form popularized by Box and Jenkins (1976). Their discussion remains one of the more complete on the subject.
To illustrate the effect of prewhitening, I will use both bidirectional and unidirectional examples. For the bidirectional example, I examine the lead/lag relationship between sea surface temperature and sea level pressure in the North Pacific. The time series I use are the first principal components for each variable derived from a PCA on monthly data for the years 1950-1992 (i.e., 516 consecutive months). These PCs describe the time trajectory of the dominant spatial modes of variability (see Figures 2.11 and 2.16 for the EOF maps). I term these "nonseasonal" PCs, as opposed to seasonal or monthly PCs which have only one value per annum.
The original time series, their autocorrelation functions, and the filtered series are illustrated in Figure 1.10. The SST time series has a particularly strong autocorrelation structure (lag 1 coefficient ~ 0.84) with significant correlation between values as far as 16 months apart. The SLP lag 1 coefficient was approximately 0.22. Both series were fit to simple AR(1) models, optimally derived via the fitting methodology described earlier:
where at is a random error term and the values in parentheses are standard error estimates of the AR coefficients.
The change in the CCF between the filtered and unfiltered time series is dramatic (Figure 1.11). Prior to prewhitening, the CCF indicated a strong bidirectional link between SLP and SST. The sign of the correlations is negative, thus anomalously low SLP results in anomalously high SST. SLP appears to lead SST at all lags between 0 and 16 months, while SST appears to lead SLP for up to 7 months. Such a relationship would be indicative of very substantial feedback and lengthy memory. The prewhitened series have a much different relationship. The prewhitened series indicate that SLP has only a lag 0 and lag 1 relationship with SST. The lag 0 cross correlation coefficient was nearly unchanged by the filtering operation (-0.350 before vs. -0.376 after). The lag 1 coefficient decreased from -0.408 to -0.226. The interpretation is that SLP, presumably through associated winds, drives variability in SST, but not vice versa. The appearance of a contemporaneous (lag 0) relationship is likely due to the lack of resolution at finer than 1 month samples. This result is very consistent with the work of Davis (1976), the first oceanographer to demonstrate that the atmosphere tends to drive the ocean and that the oceanic memory of atmospheric events is on the order of 1 month (see discussion in Frankignoul (1985)).
For the unidirectional example, I will use the same spring SST series that has been developed throughout this chapter and a catch time series for western Alaska sockeye salmon. Several authors have hypothesized that sea surface temperature may play an important role in salmon production, and this topic is further analyzed in Chapter 5. The salmon time series is described in Chapter 3, for this example it is limited to the 1950-1992 time period so as to match the SST time series.
The prewhitening model for spring SST was shown earlier to be a simple AR(1) model:
The high level of persistence in spring SST can be anticipated by the lag 12 autocorrelation coefficient for the annual SST PC. This filter was then applied to the sockeye salmon time series and the two resultant residual time series cross correlated. The before and after CCFs are given in Figure 1.12.
Prior to prewhitening, the CCF indicates a rather unusual relationship between the two series, i.e., one in which there appears to be feedback. Sea surface temperature appears to be related to salmon catch at all lags between -1 and +6 years. Following prewhitening, the number of significant lags decreases to two: lags 0 and 2 and these barely exceed the two standard error bars. If these correlations are not spurious, they indicate that western Alaska salmon production is positively related to sea surface temperature both two years prior to catch and during the year of catch. Potential mechanisms for these lag relationships are discussed in Chapter 4. I will note here that further prewhitening to compensate for a "regime effect" eliminates the significance of these lags as well.
The effects of autocorrelation on the CCF, as described in the preceding sections, have long been known, though often underappreciated. Interventions within a time series can have the same effect due to the similarity between a step intervention and large lag 1 autocorrelation coefficient. In this section, I use simulation to demonstrate how time series that have experienced an intervention effect - contemporaneous or not - will show spurious evidence of being significantly related series.
For the first simulation I generated two sets (set A and set B) of 100 time series of random deviates. Set A time series are referenced as Xt, set B time series as Yt. Each time series is of length 100, zero mean, and unit variance. In the second half of each time series, an intervention effect was added, thus the time series were of the form:
The value of the intervention effect ranged from 0.0 (no intervention effect) to 2.5 in increments of 0.5. The CCF was calculated for 100 pairs of time series (set A time series 1 with set B time series 1, and so on). At each lag, the average CCF coefficient was computed as well as the standard error (SE) of the mean. Coefficients were deemed significantly different than zero if the mean value was more than ±2*SE values from zero. The resultant CCFs for each of the 6 levels of intervention effect are illustrated in Figure 1.13. At I = 0.0 (no intervention effect), the average CCF coefficient at all lags between -6 and 6 is almost identically zero. As the intervention effect is increased, the entire set of CCF coefficients gets larger and, at I = 1.5 the entire set is significantly greater than 0. At an intervention effect of 2.5, the set of CCF coefficients ranges from 0.5-0.7.
To illustrate how large the intervention effects are, I list in Table 1.4, the average mean and variance for the 100 set B (the Yt's) time series at the different intervention levels. As expected, with no intervention effect, the average mean is zero and average variance is 1.00. The variance is doubled (i.e., the intervention effect accounts for half the total series variance) at a level of I=2.00. At this level, all CCF coefficients are highly significant, ranging around 0.4 in value.
For the second simulation, I used two sets of time series, but made the second set a lag 1 function (w1=0.5) of the first set plus additive noise, and used the same intervention rules:
The results are plotted in Figure 1.14. With no intervention
effect, only the lag 1 coefficient is statistically significant (
).
As the magnitude of the intervention effect increases, so do all the CCF
coefficients. At an intervention level of 1.5, the 2*SE bars overlap between
lag 1 and the other lag coefficients, all of which are significantly different
than zero. The average mean and variances for the 100 Yt time
series are given in Table 1.3. Once again, an intervention effect of 2.0
results in a doubling of the variance.
I experimented changing the timing of the intervention effect so that it did not coincide between the two sets of series. There was very little change in the results if the difference in timing was less than 10 time steps. This is probably due to the relatively long nature of these time series: even with a difference of 10 time steps, 80 pairs of values are matched by regime leading to the high CCF coefficients at negative and positive lags.
In this chapter, I have attempted to provide a rationale for the use of modern statistical methods in the analysis of environmental and biological data. The two principal methods I have discussed - principal component analysis and time series analysis represent a set of tools that are well suited to take advantage of the "problems" inherent in time and spatially based data. The literature is replete is with examples of physical/biological correlations that do not "stand the test of time", i.e., they are significant at the time of computation but not for long after. I surmise that these breakdowns are not due to lack of real mechanisms but rather due to misapplication of statistical procedures, either through misuse or misspecification. The single greatest problem arises due to autocorrelation and its resultant impact on the cross correlation function between series. Spatially, selection of point locations to represent larger scale variability can also be quite problematic. The use of EOF maps can assist in determining the spatial extent and regional intensity of large-scale processes.
It is my premise that fisheries-oceanography, or more generally fisheries-environmental, interactions are both real and important. In many cases, they cannot be studies via traditional techniques. That does not, however, lead directly to the dismissal of such studies. Rather, it becomes imperative that formal methods of analysis be utilized and their shortcomings acknowledged. Only in this manner will we make real progress in attempting to understand how climate and fisheries interact at the large temporal and spatial scales so important to management.
Table 1.1. The six modes by which a principal components analysis (PCA) can be conducted on a multivariate dataset. This dissertation consists solely of S-mode analyses.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Table 1.2 Univariate (U) ARIMA and Intervention (IV) models fit to dominant spring sea surface temperature principal component time series. Values in parentheses are standard errors of parameter estimates. An underlined coefficient indicates it is not significantly different than zero.
|
|
|
|
|
U1 |
(1,0,0) with constant |
SSTt = -0.011 + 0.500SSTt-1 + at (0.258) (0.136) |
|
|
|
(0.134) |
|
|
|
(0.197) (0.154) |
|
|
|
(0.152) |
|
|
|
(0.382) (0.139) (0.231) |
|
|
|
(0.141) (0.232) |
|
IV1 |
(1,0,0) no constant |
SSTt = 0.325SSTt-1 + at + 0.846It1977 (0.148) (0.296) |
|
|
|
(0.178) (0.156) (0.289) |
|
|
|
(0.150) (0.246) |
Table 1.3 Summary statistics for Univariate (U) ARIMA and Intervention (IV) models fit to dominant spring sea surface temperature principal component time series. Model abbreviations taken from Table 1.4. MAE = mean absolute error of fitted values, RSS = residual sum of squares, m = number of estimated parameters, s2a = unbiased residual variance, r2 = coefficient of determination, LL = log likelihood, AIC = Akaike's Information Criterion, SBC = Schwarz's Bayesian Criterion, and Q = portmanteau residual autocorrelation test (up to lag 10) and associated p-value.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Table 1.4 Resultant average mean and average
variance for the 100 set B simulated time series of length 100 (~N(0,1))
at different levels of an additive intervention effect.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Table 1.5 Same as Table 1.2, but the set
B time series are related to the set A time series by the function: yt
= 0.5xt-1.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|