FRACTAL ANALYSIS OF FISHERIES AND ENVIRONMENTAL TIME SERIES FOR THE DEVELOPMENT OF EARLY WARNING INDICATORS by Rodrigo Marco Montes-Aste B.Sc. University of Concepción, Chile, 1997 Marine Biologist, University of Concepción, Chile, 1997 M.Sc., University of Concepción, Chile, 2004 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Oceanography) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) October 2015 © Rodrigo Marco Montes-Aste, 2015 ii Abstract Fractal theory has been used in aquatic sciences to quantify scale-invariant relationships under the form of scaling or power laws, which represent patterns that help to understand the complex structure of exploited marine ecosystems and their populations through time. Fisheries time series are inherently variable and complex because they represent the interaction between population dynamics, oceanographic forcing’s and anthropogenic pressure (exploitation rates). Thus, to better understand the temporal dynamics of fisheries time series under the framework of fractal theory, smooth pink shrimp (Pandalus jordani) daily catch time series, daily sea surface temperature and wind stress time series from the west coast off Vancouver Island were analyzed. The ultimate goal of this thesis was the early detection of changes and critical thresholds levels over which a fishery move to a depleted or unhealthy state. This is a highly desirable objective because it gives fishery managers time to intervene in order to avoid a fishery collapse. I identified fractal patterns in smooth pink shrimp daily catches, sea surface temperature and wind stress time series. Those patterns were quantified through the Hurst coefficient (H), which integrates time series’ long-range dependence and intermittency into a single index. Inter-annual changes in Hurst coefficients from sea surface temperature appear to be related with major oceanographic anomalous conditions such regime shifts and El Niño events. Long-range dependence and intermittency were estimated for catch time series using the fractionally differenced parameter d and the Lévy index α and then combined into a single index (1/LHE indicator). Our proposed indicator, calculated using rolling windows across the history of the fishery, tracks well the vulnerable biomass for all years, which means that it is possible to have an indicator of major changes in the biomass of Pandalus jordani at the beginning and during the progress of each fishing season. iii Preface This dissertation is an original work product of the author, R.M. Montes. A version of Chapter 2 was published as: Montes, R.M, Perry, R.I., Pakhomov, E.A., Edwards, A.M. and J.A. Boutillier. 2012. Canadian Journal of Fisheries and Aquatic Sciences 69: 398-413. iv Table of Contents Abstract…………………………………………………………………………………………... ii Preface……………………………………………………………………………………………. iii Table of Contents………………………………………………………………………………… iv List of Tables……………………………………………………………………………………... viii List of Figures……………………………………………………………………………………. ix Acknowledgements…………………………………………………………….………………… xv Dedication………………………………………………………………………………………… xvi CHAPTER 1 Introduction to fractals and statistical considerations……………………... 1 1.1 Definitions and description of deterministic fractals……………………………………… 1 1.2 Estimation of fractal parameters…………………………………………………………… 3 1.3 Fractals in the time domain………………………………………………………………... 5 1.3.1 Fractal functions.............................................................................................................. 5 1.3.2 Long-range persistence……………………………………………..…………………. 6 1.4 Multifractals……………………………………………………………..……………….... 9 1.4.1 Definition and description of multifractals...................................................................... 9 1.4.2 Multifractals in the time domain………………………………………………...…….. 9 1.4.3 The Universal Multifractal Model (UMM)…………………………………………..... 11 1.5 Fractal analysis in marine ecology……………………………………………………..…. 12 1.6 Spatial fractal analysis in fisheries oceanography..……………………………………..… 14 1.7 Temporal fractal analysis in fisheries oceanography…………………………………….... 16 1.7.1 Case study 1 (Ferriere and Cazelles 1999)…...………………………………………... 17 1.7.2 Case study 2 (Montes 2004)…………………............................................................... 17 1.7.3 Case study 3 (Halley and Stergiou 2005)……............................................................... 19 1.7.4 Case study 4 (Niwa 2006, 2007)…..…………………………………………………... 20 1.7.5 Case study 5 (Mendes et al. 2014)……………………………………………………... 21 1.8 General objectives………………...…………………………………………………….… 22 v CHAPTER 2 Multifractal patterns in the daily catch time series of smooth pink shrimp (Pandalus jordani) from the west coast of Vancouver Island, Canada…………………………………………………… Vancouver Island, Canada ………………………........................................... 23 2.1 Introduction…………………………………………………………………………….….. 23 2.2 Methods…………………………………………………………………………….…….... 26 2.2.1 Fractals and long-memory processes…………………………………………....…….. 26 2.2.2 Pink shrimp fishery and total daily catches………………………………………...….. 31 2.2.3 Monofractal and multifractal analysis…….………………………………………….... 37 2.2.4 Estimation of confidence intervals…….......................................................................... 41 2.2.5 Detection of scaling ranges……………………………………………………………. 41 2.2.6 Selection between monofractal and multifractal models……………………………… 42 2.3 Results……………………………………………………………………………………... 43 2.3.1 Scaling TDC and stationarity of TDCr time series……………………………………. 43 2.3.2 Selection of scaling ranges…………………………………………………………...... 46 2.3.3 Estimation of monofractal and multifractal parameters and selection of models………………………………………………………..………........ 48 2.4 Discussion…………………………………………………………………………………. 50 2.4.1 Identification of multifractal patterns…………………………………………………. 50 2.4.2 Potential applications of fractal models in the pink shrimp fishery ………………...… 52 2.4.3 Hypotheses for the causes of these fractal patterns……………………………………. 55 CHAPTER 3 Estimation of long-range temporal dependence of sea surface temperature and wind stress and their effects on smooth pink shrimp (Pandalus jordani) fisheries catches off the west coast of Vancouver Island………………………………………….. 60 3.1 Introduction……………………………………………………………………………….. 60 3.2 Methods…………………………………………………………………………………… 63 3.2.1 Oceanographic databases……………………………………………………………… 63 3.2.2 Fishery database………………………………………………………………………. 66 3.2.3 Removal of unwanted oscillations…………………………………………………….. 72 3.2.4 Fractal analysis ………………………………………………………………….…….. 73 vi 3.2.5 Multifractal analysis…………………………………………………………………… 76 3.2.6 Bivariate multifractal analysis………………………………...…………………….… 77 3.3 Results…………………………………………………………………………………….. 79 3.3.1 Filtered sea surface temperature and wind stress time series…………………………. 79 3.3.2 Long-range dependence of oceanographic and fishery time series…………………… 81 3.3.3 Multifractality in oceanographic and fishery time series……………………………… 87 3.3.4 Bivariate long-range dependence between oceanographic and fishery time series…… 87 3.4 Discussion…………………………………………………………………………………. 89 3.4.1 Long-range dependence in oceanographic time series………………………………… 89 3.4.2 Variability and representativeness of Hurst coefficients……………………………… 91 3.4.3 Selection of a long-range dependence over a short-range dependence model………… 93 CHAPTER 4 A self-similar early warning indicator accounting for long-range dependence and intermittent variability in Pandalus jordani fishery catches from the west coast of Vancouver Island……………..… 97 4.1 Introduction…………………………………………………………………………...…… 97 4.2 Methods……………………………………………………………………………………. 100 4.2.1 Estimation of Hurst coefficient confidence intervals for simulated gappy data….…… 101 4.2.2 Complexity of TDC time series and removal of unwanted oscillations………….…… 103 4.2.3 Estimation of Hurst coefficient, LRD and intermittency parameters………………….. 112 4.2.4 Comparison between LHE’s and classical statistical indicators……………………… 114 4.3 Results…………………………………………………………………………………….. 116 4.3.1 Hurst coefficient confidence intervals for simulated gappy data…………………….. 116 4.3.2 Local Hölder exponents, LRD and intermittency parameters………………………… 122 4.3.2.1 Early detection of changes in total daily catches (TDC)………………………….. 122 4.3.2.2 Early detection of changes in vulnerable biomass………………………………… 126 4.3.3 Comparison of early warning indicators…………………………………..………….. 130 4.4 Discussion…………………………………………………………………………………. 133 4.4.1 Early detection of fisheries changes…………………………………………………... 133 4.4.2 Long-range correlations in fish recruitment and local Hölder exponents (LHE’s)…… 136 4.4.3 Comparison of early warning indicators………………………………………………. 138 vii CHAPTER 5 Summary and Conclusions…………………………….………………….. 142 REFERENCES ……………………………………….………………………………………. 150 APPENDIX A ………………………………………………………………………………. 190 A.1 Probability distribution for filtered sea surface temperature time series…………..…….. 190 A.2 Probability distribution for filtered wind stress time series……………………………… 191 A.3 Maximum likelihood estimation of Hurst coefficients for filtered SST time series………. for filtered SST time series 192 A.4 Filtered sea surface temperature for the 1940-2013 period……...……………………….. 193 A.5 Diagram showing stages and decisions adopted for the estimation of Hurst and generalized Hurst coefficients in univariate and bivariate fractal analysis….…………..… 194 A.6 Probabilities of rejecting the null hypothesis when using a long-range dependent (LRD) model to describe a short-range dependent (AR1) process……………………………….. 195 viii List of Tables Table 2.1 Definitions of variables and parameters used in Chapter 2…………………... 28 Table 3.1 Definitions of variables and parameters used in Chapter 3……………….….. 68 Table 4.1 Definitions of variables and parameters used in Chapter 4…………………... 105 ix List of Figures Figure 1.1 Simulated fractional Gaussian noises (fgn; left column) and fractional Brownian motions (fBm; right column) for increasing levels of Hurst coefficients H=0.5, 0.6, 0.7, 0.8, 0.9 for time series of length N=500. Simulations were conducted in R software (R Development Core Team 2005, version 3.0.2) according to the method described in Davies and Harte (1987) using the R package “fracdiff” (Fraley et al. 2006)…………….……………... 8 Figure 2.1 Fisheries and Oceans Canada Pacific Fishery Management Areas (PFMA). Data for this study came from PFMA 121, 123, 124 and 125 off the west coast of Vancouver Island, Canada…………………………………………… 32 Figure 2.2 Pink shrimp (Pandalus jordani) total daily catch (i.e. daily sum of all recorded catches; TDC) time series from 1994 to 1996 for the April-October period (642 fishing days) in Pacific Fishery Management Areas 121, 123, 124, 125 west of Vancouver Island, Canada. A sinusoidal curve with an annual period of 214 fishing days was fitted to this time series. A slightly increasing linear trend is also visible in this time series………………………. 33 Figure 2.3 Diagram showing the transformations and mathematical operations conducted on pink shrimp (Pandalus jordani) total daily catch time series (TDC) on a step by step basis. TDCr represents the total daily catch residual time series and TDCa is the total daily catch anomaly time series obtained after cumulatively summing the TDCr time series……………………………. 36 Figure 2.4 Power spectrum of the pink shrimp (Pandalus jordani) total daily catch time series (TDC) off the west coast of Vancouver Island for 1994-1996. The fitted straight line is a preliminary indication of scaling……………………… 44 Figure 2.5 Structure function Yr(q) versus r in a log-log plot obtained from the pink shrimp (Pandalus jordani) total daily catch residual time series (TDCr) for the west coast of Vancouver Island from 1994-1996. Two increasing trends, the first between 1 and 3 fishing days and the second between 4 and fishing 15 days are visible…………………………………………………………….. 45 Figure 2.6 Scaling range (between 16 and fishing 120 days) selected by the adaptive search algorithm (Xia et al. 2003) applied to detect the linear region of the structure function curve of the total daily catch anomaly time series (TDCa) for q=1. Upper scaling limit is marked by an arrow. Inset shows a close-up view of the departure of the structure function from the linear trend at 120 fishing days……………………………………………………………………. 46 Figure 2.7 Temporal variability of local slopes calculated for TDC anomaly time series (TDCa) for time lags between 16 and 321 fishing days. Local slopes fluctuate around a constant value up to approximately 120 fishing days and notably decrease after that……………………………………………………………... 47 x Figure 2.8 Structure function Yr(q) versus r in a log-log plot obtained from the pink shrimp (Pandalus jordani) total daily catch anomaly time series (TDCa) for the west coast of Vancouver Island from 1994-1996. A linear trend selected by the adaptive search algorithm (Xia et al. 2003) for q=1 in Figure 2.6 is visible also for statistical moments q=1, 2, 3, 4 for a range of scales from 16 to 120 fishing days marked between vertical dotted lines…………………….. 49 Figure 2.9 Empirical curve for scaling exponents, ζ(q), (filled circles) compared to the monofractal curve ζ(q)=qH (continuous line) and to the theoretical universal multifractal function (crosses) obtained from the pink shrimp (Pandalus jordani) total daily catch anomaly (TDCa) time series with H=0.8779, C1= 0.0429 and α=1.5073…………………………………………………………. 50 Figure 3.1 Amphitrite Point Lighthouse (48º55’N; 125º32’W, circle) and buoy number 206 (48.84 ºN; 126º00’W, triangle) located off the west coast of Vancouver Island, Canada………………………………………………………………… 64 Figure 3.2 Sea surface temperature time series temperature (SST in °C) at Amphitrite Point Lighthouse for the 1940-2013 period. Missing values (3.6%) were replaced by a long-term mean (see text) and their positions marked by red crosses centered at 5ºC. A horizontal line is centered at time series mean (10.10ºC). Vertical grey lines separate periods of ten years………………….. 65 Figure 3.3 Time series of daily mean wind stress (WS, m/s2) for the 1958-1998 period and its mean annual cycle (green line). b) Residual wind stress time series. Vertical lines separate periods of ten years…………………………………… 66 Figure 3.4 a) Time series of smooth pink shrimp total daily catches (TDC, 1000 Kg) for the 1987-1998 period, and b) TDC anomaly time series. Missing values (16.5%) were replaced by a long-term mean (see text) and their positions marked by red crosses centered at 1000 Kg. Vertical grey lines separate periods of 428 days, which correspond to two continuous fishing seasons…... 67 Figure 3.5 Sea surface temperature wavelet power spectrum for 1940-2013 period. Morlet wavelet was selected as mother wavelet. A white dotted line indicates the cone of influence (COI), below which power spectrum coefficient should be interpreted with caution. A frequency band in dark red associated to the annual cycle is clearly visible……………………………………………….... 79 Figure 3.6 a) Time series of normalized sea surface temperature (SST in °C) for the 1940-2013 period, and b) Wavelet reconstructed sea surface temperature representing variability between 128 and 4096 days (0.35 and 11.22 years, respectively). The period June 1997- June 1998 (highlighted with a orange bar) represents one of the two strongest El Niño events detected during the last century. Vertical lines separate periods of ten years……………………… 80 xi Figure 3.7 a) Time series of normalized daily mean wind stress (WS, m/s2) for the 1958-1998 period, and b) wavelet reconstructed wind stress representing variability between 128 and 4096 days (0.35 and 11.22 years, respectively). Vertical grey lines separate periods of ten years……………………………………….. 81 Figure 3.8 Hurst coefficients estimated for filtered sea surface temperature (SSTF) time series during the 1940-2013 period. The horizontal dotted line indicates the mean of Hurst coefficients (0.9661). Vertical grey lines separate periods of ten years. Major El Niño events (1957, 1983, 1998) are demarked with filled red circles. Year 2011 (filled orange circle) coincides with an alternation of strong El Niño/La Niña events. Regime shifts detected in the North-East Pacific for 1988 (Hare and Mantua 2000) and 1975 (Overland et al. 2008) were highlighted with a filled blue circle……..………………………………… 83 Figure 3.9 Root mean square fluctuation function F(s) for filtered sea surface temperature time series (SSTF). a) Scaling range (blue circles) extends between 16 and 1024 days (24-210, H=0.8760). A scaling break marked with an arrow is visible around 256 days. b) Scaling range (blue circles) extends between 16 and 80 days (24-26.32, H=0.9675). Lines were fitted to scaling ranges to estimate the Hurst coefficient……………………………………….. 84 Figure 3.10 Root mean square fluctuation function F(s) for filtered wind stress time series (WSF) for the 1958-1998 period. a) First scaling range (grey line fitted to blue circles) extends between 16 and 120 days (24-26.90) (H=0.7538). b) Second scaling range (grey line fitted to red circles) extends between 512 and 3726 days (H=0.6563)…………………………………………..…………………... 85 Figure 3.11 Root mean square fluctuation function F(s) for smooth pink shrimp total daily catch (TDC) anomaly time series for the 1987-1998 period. a) First scaling range (grey line fitted to blue circles) extends between 16 days and 120 days (24- 26.90, H=0.8568). b) Second scaling range (grey line fitted to red circles) extends between 121 and 642 days…………………………….…. 86 Figure 3.12 Bivariate long-range cross-correlations between filtered sea surface temperature and TDC anomaly time series for the 1987-1998 period. Scaling range (highlighted with filled circles) used to estimate bivariate Hurst exponent extends between 16 days and 80 days (24- 26.32, hof=0.8580)………. 88 Figure 3.13 Bivariate long-range cross-correlations between filtered wind stress and smooth pink shrimp (TDC) anomaly time series for the 1987-1998 period. Scaling range (highlighted with filled circles) used to estimate the bivariate Hurst exponent extends between 16 and 120 days (24- 26.90, hof=0.8315)…….. 89 xii Figure 4.1 Wavelet power spectrum for smooth pink shrimp total daily catch (TDC) time series for 1987-1998 fishing seasons. Vertical axis represents scales in days. Shaded area between 128 and 512 days contains most of the total daily catch time series variance. The Morlet wavelet was selected as the mother wavelet. Black contours encompass regions for which the lag-1 autoregressive coefficient (AR) is significantly different from 0.65. White dotted lines indicates the cone of influence (COI), below which the power spectrum coefficient should be interpreted cautiously………..………….……………… 108 Figure 4.2 a) Time series of smooth pink shrimp total daily catch (TDC) for 1987-1998 fishing seasons, b) catch level shifts from 128 to 512 days (step-like changes) were calculated using Haar wavelets. Total daily fishing effort (TDE in trawl hours, green line) follows major step-like changes in TDC, and c) TDC residuals obtained after subtracting step like changes in catches obtained in (b). Vertical grey lines separate annual fishing seasons………………………. 110 Figure 4.3 Histograms for the smooth pink shrimp total daily catch (TDC) fishing seasons from 1987 to 1998……………………………………………………. 111 Figure 4.4 a) Standardized smooth pink shrimp total daily catch (TDC) time series from 1987 to 1998. Vertical grey lines separate annual fishing seasons. b) Histogram of standardized TDC time series……………………………………………...….. 116 Figure 4.5 Uncertainty and bias in the estimation of H from simulated gappy time series using the quantile/maximum likelihood method (QQ-ML) and the KG method as described in the text for α=1.6 and six increasing long-range dependence levels d. Red horizontal lines mark theoretical H value. Experimental H values observed in boxplots were calculated from replicates obtained using ME bootstrap………………………………………………….. 117 Figure 4.6 Uncertainty and bias in the estimation of H from simulated gappy time series using the quantile/maximum likelihood method (QQ-ML) and the KG method as described in the text for α=1.7 and six increasing long-range dependence levels d. Red horizontal lines mark theoretical H value. Experimental H values observed in boxplots were calculated from replicates obtained using ME bootstrap………………………………………………….. 118 Figure 4.7 Uncertainty and bias in the estimation of H from simulated gappy time series using the quantile/maximum likelihood method (QQ-ML) and the KG method as described in the text for α=1.8 and six increasing long-range dependence levels d. Red horizontal lines mark theoretical H value. Experimental H values observed in boxplots were calculated from replicates obtained using ME bootstrap………………………………………………….. 119 xiii Figure 4.8 Uncertainty and bias in the estimation of H from simulated gappy time series using the quantile/maximum likelihood method (QQ-ML) and the KG method as described in the text for α=1.9 and six increasing long-range dependence levels d. Red horizontal lines mark theoretical H value. Experimental H values observed in boxplots were calculated from replicates obtained using ME bootstrap………………………………………………….. 120 Figure 4.9 Uncertainty and bias in the estimation of H from simulated gappy time series using the quantile/maximum likelihood method (QQ-ML) and the KG method as described in the text for α=2.0 and six increasing long-range dependence levels d. Red horizontal lines mark theoretical H value. Experimental H values observed in boxplots were calculated from replicates obtained using ME bootstrap………………………………………………….. 121 Figure 4.10 Variability of local Hölder exponents (LHE’s, increasing red line) and total daily catch residual time series (black lines) for 1993-1998 fishing seasons…. 123 Figure 4.11 Variability of local d’s (long-range dependence parameter, increasing red line) and total daily catch residual time series (black lines) for 1993-1998 fishing seasons………………………………………………………………… 124 Figure 4.12 Local slopes estimated for 1993-1998 fishing seasons using the method described by Takalo (1995) and Maraun et al. (2004). A window equal to 91 data points was used for the estimation of d’s (long-range dependence). Major changes in slopes are visible for the 1995 fishing season……………… 124 Figure 4.13 Variability of local α (intermittency parameter, red line) and total daily catch residual time series (black line) for 1993-1998 fishing seasons………………. 125 Figure 4.14 (Upper plot) Time series of total daily catch (TDC) residuals for 1993-1998 fishing seasons. (Lower plot) Inverse Local Hölder Exponent (1/LHE’s, red line) estimated using the KG method is presented for a better visualization of early changes in exponents in relation to biomass fluctuations. Vertical lines separate annual fishing seasons. Blue squares represent vulnerable biomass reconstructed for Pacific Fisheries Management Areas 124 and 125 off the west coast of Vancouver Island (e.g. Fig. 3.1) according to Martell et al. (2000). Shaded area encloses fishing seasons 1994 and 1995, for which major changes in biomass occurred. Biomass levels were localized on the plot at the end of April for each fishing season because fishery independent biomass surveys were conducted at the end of April or beginning of May each year…. 127 Figure 4.15 (Upper plot) Time series of total daily catch (TDC) residuals for 1993-1998 fishing seasons. (Lower plot) Inverse local d parameter (1/d, red line) with vertical grey lines separating annual fishing seasons. Blue squares represent vulnerable biomass reconstructed for Pacific Fisheries Management Areas 124 and 125 according to Martell et al. (2000). Shaded area encloses fishing seasons 1994 and 1995, for which major changes in biomass occurred. Biomass levels were localized on the plot at the end of April for each fishing season because fishery independent biomass surveys were conducted at the end of April or beginning of May each year…………………………………... 128 xiv Figure 4.16 (Upper plot) Time series of total daily catch (TDC) residuals for 1993-1998 fishing seasons. (Lower plot) Local intermittency exponent (α, red line) with vertical grey lines separating annual fishing seasons. Blue squares represent vulnerable biomass reconstructed for Pacific Fisheries Management Areas 124 and 125 according to Martell et al. (2000). Shaded area encloses fishing seasons 1994 and 1995, for which major changes in biomass occurs. Biomass levels were localized on the plot at the end of April for each fishing season because fishery independent biomass surveys were conducted at the end of April or beginning of May each year………………………………………….. 130 Figure 4.17 Dynamics of local Hölder exponents (LHE, red line) and fluctuations of total daily catch (TDC) residuals for 1993-1998 fishing seasons (black continuous line) in comparison to a) autoregressive lag-1 coefficient, b) density ratio fluctuations, c) skewness and d) kurtosis (blue lines). Vertical grey lines separate annual fishing seasons……………………………………………….. 131 Figure 4.18 Patterns of change of local slopes for local Hölder exponents (LHE’s, continuous red line) in comparison to a) local slopes of autoregressive lag-1 coefficients (continuous grey line) and local slopes of density ratios (dotted grey line); and b) local slopes of skewness (continuous grey line) and local slopes of kurtosis (dotted grey line) (scaled to fit this figure). Vertical grey lines separate annual fishing seasons. A window equal to 91 observations was used to calculate slopes……….……………………………………………….. 132 Figure 4.19 (Upper plot) Inverse local Hölder exponents (1/LHE, red line) and vulnerable biomass for PFMA 124 and 125 (Martell et al. 2000). (Middle plot) Recruitment anomalies estimated for PFMA 124. (Lower plot) Recruitment anomalies estimated for PFMA 125 (Martell 2002)…………………………... 138 xv Acknowledgements This thesis represents the culmination of an idea that starts many years ago before my arrival to The University of British Columbia, Vancouver. My supervisors at the Department of Earth, Ocean and Atmospheric Sciences, Dr. Evgeny Pakhomov, and my co-supervisor at the Department of Fisheries and Oceans (DFO) Canada, Dr. Ian Perry, have significantly contributed to transform this idea into this tangible work. Both of them strongly supported me, especially during the last period of time when life became difficult for me and for all my family members. I really owe you my most sincere “THANKS” for your invaluable help. I also express my gratitude to my committee members, Dr. William Hsieh from Earth, Ocean and Atmospheric Sciences (UBC) and Dr. Steve Martell from Fisheries Centre (UBC) for their invaluable help. Dr. Andrew Edwards, from Department of Fisheries and Oceans (DFO) Canada, was always present to give me an advice to improve this work. I am really grateful to the Department of Earth, Ocean and Atmospheric Sciences and to the Faculty of Graduate Studies from The University of British Columbia who allowed me to finish my thesis work, and consequently, to go one step ahead in my scientific career. My friend, Dr. Renato Quiñones supported me, especially during stressful moments. His clear ideas on the resolution of theoretical concepts contributed to finishing this thesis. Last but not least, I also note the work from the staff of The University of British Columbia Student Housing and Hospitality Services who made my stay in Acadia Park, along with my family members, a beautiful and incredible dream. xvi Dedication I would like to dedicate this piece of work to my daughters Pia, Isidora and Martina who are an essential part of my family, and to my parents, Marco and Eliana, who always encourage me to do what I really want to do; no matter how difficult or hard it looks. My sister Carla and my brother Javier were always present and accompanying me especially when I was going through a rough time. 1 Chapter 1: Introduction to fractals and statistical considerations 1.1 Definition and description of deterministic fractals Geometrical objects can be classified according to their characteristic length, which is the typical length of an object (e.g. the radius of a sphere, the diameter of a circle, the height of a tower). All objects of Euclidean geometry possess their own characteristic lengths and the shapes of such objects are mostly smooth, i.e., in the majority of their surface points a derivative can be determined. In other words, their surface is differentiable almost everywhere (Ficker & Benešovský 2002). Nevertheless, there are other objects besides Euclidean objects, which are not differentiable in any of their surface points because their structure lacks a characteristic length. With these objects it is no longer possible to find a length scale, which would show characteristic structural features because the same geometrical structures are present on all their length scales (Mandelbrot 1983). If the observational scale, under which we are analyzing these objects, is reduced, new components arise, and the new ones are similar to those of the original unit (Evertsz and Mandelbrot 1992). These special classes of objects are called fractals, and in this context a fractal can be defined as an object composed of sub-units (and further sub-units) that resemble the larger scale structure, a property known as self-similarity (Mandelbrot 1983, Feder 1988, Goldberger et al. 2002). This property expresses that each part of a fractal object is a linear geometric reduction of the whole, with the same reduction ratios in all directions (in contrast to self-affine processes in which the reductions are still linear but the reduction ratios are different in different directions). Self-similar processes are invariant in distribution under scaling of time and/or space, and the scaling coefficient (or index of self- 2 similarity) is a non-negative number denoted H or the Hurst coefficient1 (Samorodnitsky & Taqqu 1994). Due to this property, the shape of fractals is nonrectifiable2, consisting of an infinite sequence of clusters within clusters or waves within waves3. In rectifiable objects, increasingly accurate measurements based upon successive scale reductions converge to a limit, which is the true extent (length) of the object (Li 2000). However, in fractals4, the same procedure generates infinite series according to a power law: 1.1 𝜆 𝑢 = 𝑘 𝑢!!! here λ is the length of the object measured at the scale unit u (the length of the object diverges as u → 0), k is a constant, and D is the fractal dimension. This latter dimension allows one to describe conceptual or concrete objects that realize "a certain degree" of occupation of a bi- or tri-dimensional Euclidean space, somewhere between a point and a line, between a curve and a surface, or between a surface and a volume. D has to be considered as a measure of that degree of occupation (Frontier 1987), and increases with the increase in complexity of the geometric object (Mandelbrot 1977, 1983). The fractal dimension can be estimated after taking logarithms5 of both sides of equation 1.1, so it can be expressed as: 1.2 log 𝜆 𝑢 = log 𝑘 + 1− 𝐷 log 𝑢 Thus, D can be estimated by fitting a linear regression between the measurements of the 1A process 𝑋 = {𝑋 𝑡 , 𝑡 ∈ ℝ} is self-similar with index H if, for any a>0, the finite dimensional distributions of {𝑋 𝑎𝑡 , 𝑡 ∈ ℝ} are the same as those of {𝑎!𝑋 𝑡 , 𝑡 ∈ ℝ}. 2 The term nonrectifiable means that they are infinite in length (Steinhaus 1954). 3 In other words, a fractal outline is a jagged, non-differentiable curve, often associated with a self-similar structure that is repeated at different scales (Hastings and Sugihara 1993). 4 Stochastic fractals (in contrast to deterministic fractals) are self-similar only in a statistical sense (Voss 1985). 5 Logarithm to the base e (ln), logarithm to the base 10 (log) and logarithm to the base 2 (log2) are indistinctively used in fractal research. In almost all applications we will require the ratio of logarithms, and in consequence the results are the same if for example log or ln is used (Turcotte 1997). 3 different lengths (λ) obtained through varying the scale unit u. The fractal concept can be applied not only to irregular geometric forms that lack a characteristic (or single) scale of length, but also to complex processes that lack a single time scale. Fractal (or scale-invariant) processes generate irregular fluctuations6 on multiple time scales, analogous to fractal objects that have a wrinkle structure on different length scales (Goldberger 1996). In an idealized model, scale-invariance holds on all scales, but the real world imposes upper and lower bounds (smax and smin, respectively) over which such behavior applies (Goldberger et al. 2002). For these types of fractals the scaling range (SR) is given in decades as: 1.3 𝑆𝑅 = 𝑙𝑜𝑔!"(𝑠!"# 𝑠!"#) Some authors argue that the minimum scaling range needed in order to affirm that the analyzed signal is fractal is two decades (e.g., Avnir et al. 1998; Biham et al. 1998a,b; Mandelbrot 1998) but this is a controversial issue that is still under discussion. 1.2 Estimation of fractal parameters The estimation of fractal parameters should be conducted considering statistical properties of the time series (stationarity vs. non-stationarity) and the category of the method (parametric, non-parametric, semi-parametric and heuristic) (Beran et al. 2013). It is important to remark that, in contrast to deterministic fractals presented in section 1.1, this thesis is focused on statistical fractals. This means that the fractal properties detected in analyzed time series are based on an ensemble average, which is the statistical average of multiple time series realizations. For the ensemble average, each realization is expected to 6 Although fractals are irregular, not all irregular structures or time series are fractal. 4 break the scaling in some way and also to recover only if sufficient data points are averaged (Lewis et al. 2004). As we usually work with only one time series, long time series are required to obtain reliable estimators of fractal parameters (Cannon et al. 1997). In consequence, we are assuming that time series are ergodic, which means that time averages for a unique time series realization converge to the ensemble average (Chatfield 2000). Many statistical methods are available in the scientific literature for the estimation of the fractal dimension (D), the Hurst parameter (H) and the long-range dependence parameter (d); the latter being directly related to H (Samorodnitsky and Taqqu 1994). Most recognized methods were developed for the estimation of fractal parameters assuming stationarity; i.e. the mean and variance of the time series do not vary substantially over time. Into this category, the most recognized parametric and semi-parametric methods are Gaussian maximum likelihood and Whittle estimation methods (Pilgram and Kaplan 1998, Faÿ et al. 2009, Robbertse and Lombard 2012). The non-parametric methods include Log-Periodogram Regressions, Local Whittle and Log-Wavelets Regression methods (Beran et al. 2013), among others. On the other hand, the heuristic methods developed for the analysis of stationary time series are numerous, being Dispersional Analysis (Caccia et al. 1997, Cannon et al. 1997), Semi-variograms (Burrough 1981, Dolan et al. 1998), Structure Functions (Tennekoon et al. 2005, Lovejoy and Schertzer 2012), Power Spectrum (Feder 1988, Witt and Malamud 2013), Scaled Windowed Variance (Mandelbrot 1985, Peng et al. 1994), Rescaled-Range Method (Mandelbrot and Wallis 1969, Davies and Harte 1987) and Detrended Fluctuation Analysis (DFA, Peng et al. 1994, Bunde et al. 2000, Kantelhardt et al. 2001) among the most used. All of the methods used in this thesis can handle the non-stationary features of 5 time series (e.g. DFA, Chapter 3) or a filter was applied for the elimination of trends and periodicities (e.g. wavelets, Chapters 3 and 4) or any other component that makes a time series non-stationary. 1.3 Fractals in the time domain In this section we briefly define and explain the two most common stochastic fractal functions and their basic mathematical properties. 1.3.1 Fractal functions The two basic types of stochastic fractal functions are characterized by their long-term memory content and self-similarity structure: (i) fGn (fractional Gaussian noise) and (ii) fBm (fractional Brownian motion) (Taqqu et al. 1995). The fGn time series corresponds to a stationary stochastic process consisting of a set of normally distributed random variables with zero mean where correlation between nth neighbors separated in time by τ = nΔt is given by: 1.4 𝜌! = 12 |𝜏 + 1|!! − 2|𝜏|!! + |𝜏 − 1|!!) where 𝜌! is the correlation coefficient between observations separated by a lag time equal to τ (Caccia et al. 1997). As we previously explained, the Hurst coefficient (H) is a measure of the long-range dependence of a self-similar series (e.g., Taqqu et al. 1995). When H> 0.5, all the observations of a time series classified as fGn are positively correlated, and the closer H is to 1 the smoother the function. When H<0.5, all observations of the analyzed fGn time series separated by a τ≥1 are negatively correlated. A time series having the characteristics of Gaussian white noise corresponds to an fGn series with H=0.5 and 6 𝜌!= 0 for τ≥1 (Sugihara and May 1990). The second fractal time series, closely related to the first, is fBm, a non-stationary time series with stationary increments (Cannon et al., 1997); this is a special case and belongs to a general family of scaling models known as fLm (fractional Lévy motion) (Taqqu 1987). The fBm time series is defined as a continuous function m(x) with the following statistical properties: 1.5 𝐸{ 𝑚(𝑥)!! −𝑚(𝑥)!!] = 0 1.6 𝐸{ 𝑚(𝑥)!! −𝑚(𝑥)!!]! = 𝑘 𝑚(𝑥)!! −𝑚(𝑥)𝑡!|!! where E{m(x)} is the mean or expected value of the temporal variable x, k is a constant of proportionality, and H is the Hurst coefficient, described earlier. The expected value is obtained by averaging over many increments of size |t2-t1|. The latter equation guarantees a fractal nature for fBm, and since it holds only on the average, statistical fractals are stochastic (Voss 1985). For the special case when H=0.5, equation (1.6) becomes: 1.7 𝐸{ 𝑚(𝑥)!! −𝑚(𝑥)!!]! = 𝑘 𝑚(𝑥)!! −𝑚 𝑥 !!| used for defining Brownian movement (cBm) in classic physics (Mandelbrot 1983). This indicates that, lacking a correlation between observations separated by different time lags, any value adopted by the time series in the future is equally likely (Sugihara and May 1990). 1.3.2 Long-range persistence Long-range persistence, also called long-range dependence, long-memory or long-range correlations, is an expression of self-similar patterns in the time domain. In this context, the extension of the idea of self-similarity in spatial objects to time series was conducted by 7 Mandelbrot and van Ness (1968) which called them self-affine fractal or self-affine time series when after appropriate rescaling of axes a statistically similar time series is produced (Witt and Malamud 2013). In this text, the term long-range dependence (LRD) was used as an expression of time series self-similarity. Is important to keep in mind that the strength of long-range dependence is directly related to the fractal dimension D and the Hurst coefficient H. A long-range dependence time series can be expressed as (Beran et al. 2013): 1.8 𝐶(𝜏) ~𝜏!(!!!), 𝜏 → ∞,−1 < 𝛽 < 1 where C corresponds to the autocorrelation function measured at lag τ. The strength of long-range dependence is estimated by the parameter β, where β=0 corresponds to a time series with no long-range dependence, β>0 to a long-range dependence time series and β< 0 to a long-range antipersistent time series. For Gaussian stationary time series (fractional Gaussian noises, -1<β<1) the Hurst coefficient 𝐻 = 𝛽 + 1 2. Simulated fractional Gaussian noises (fgn) and fractional Brownian motions (fBm) of length N=500 for increasing levels of Hurst coefficients H=0.5, 0.6, 0.7, 0.8, 0.9 are shown on figure 1.1. The first and last rows represent time series with no long-range dependence (β=0, H=0.5) and pronounced long-range dependence (β=0.4, H=0.9), respectively. 8 Figure 1.1 Simulated fractional Gaussian noises (fgn; left column) and fractional Brownian motions (fBm; right column) for increasing levels of Hurst coefficients H=0.5, 0.6, 0.7, 0.8, 0.9 for time series of length N=500. Simulations were conducted in R software (R Development Core Team 2005, version 3.0.2) according to the method described in Davies and Harte (1987) using the R package “fracdiff” (Fraley et al. 2006). 9 1.4 Multifractals 1.4.1 Definition and description of multifractals In contrast to monofractals (or simple fractals) which are homogeneous in the sense that they have the same scaling properties characterized by only one exponent or coefficient (i.e., H) throughout the entire signal, multifractals are a class of signals which requires a larger, and theoretically infinite, number of indices to characterize their scaling properties (Harte 2001). Multifractal signals are intrinsically more complex and inhomogeneous than monofractals, i.e., different parts of the signal may have different scaling properties (Goldberger et al. 2002). Whereas fractal analysis looks at the geometry of a pattern (e.g., presence/absence, true/false), multifractal analysis looks at the arrangements of quantities (e.g., population density), or in other words, a fractal can be defined as a geometrical set of points and a multifractal as a mathematical measure (Lavallée et al. 1993). The multifractal approach implies that a statistically self-similar measure can be represented as a combination of interwoven fractal sets with corresponding scaling exponents7. A combination of all the fractal sets produces a multifractal spectrum that characterizes variability and heterogeneity of the analyzed variable (Kravchenko et al. 1999). 1.4.2 Multifractals in the time domain Analyzing a multifractal process φλ under a scaling ratio of λ, we can establish that: 7 Another explanation of the relation between a fractal and a multifractal was given by Halley et al. (2004) who state that: "multifractal analysis describes an object as the sum of many fractal subsets; the sets of points which each particular value of density forms its own fractal pattern". 10 1.9 Pr 𝜙 ≥ 𝜆! ≈ 𝜆!! ! In this equation, γ is a variable known as the order of singularity, which represents a defined intermittent level (Seuront et al. 1999). The term c(γ), the scaling exponent of the probability distribution in equation (1.9), correspond to the co-dimension of the multifractal process. c(γ) represents a non-linear and concave function of γ and is given by: 1.10 𝑐 γ = 𝑑 − 𝐷(𝛾) where D(γ) characterizes the hierarchy of fractal dimensions8 associated with different intermittency levels (i.e., singularities) and d is the dimension of the space considered (d = 1 for time series) (Schertzer and Lovejoy 1987, Seuront et al. 1999). For a given scaling ratio (λ), the probability (Pr) of finding values belonging to the multifractal process that exceed an established threshold (i.e., λγ) decreases as γ increases. Thus, the probability distribution of the field φλ is scaled based only on the scale ratio (i.e., there is no characteristic scale; Boufadel et al. 2000). When c(γ) varies as a function of γ, the exceedance probability at each threshold value scales differently than another value; this property is known as multiple scaling or multiscaling. Multifractals, therefore, are characterized by an infinite hierarchy of subsets, each having its own characteristic fractal dimension (Feder 1988). In general, each subset of the hierarchy corresponds to a fraction of the space made up of observations that surpass the given threshold (Seuront et al. 1999). In the case of fractals (or monofractals), c(γ) is constant, producing a single scaling or monoscaling (Boufadel et al. 2000). Instead of characterizing the system’s statistical 8 This definition is equivalent to those given above which refers to multifractals as the "sum of many fractal subsets…" or as a "combination of interwoven fractal sets with corresponding scaling exponents". 11 properties through the scaling of its probabilities c(γ), we can do an equivalent analysis of the scaling of its statistical moments (Tessier et al. 1993, Boufadel et al. 2000). 1.4.3 The Universal Multifractal Model (UMM) There are various mathematical models for characterizing and quantifying a multifractal process (e.g. the p-model, Meneveau and Sreenivasan 1987; the log-Poisson model, She and Waymire 1995) but a general one is Schertzer and Lovejoy’s Universal Multifractal Model (UMM, Schertzer and Lovejoy 1987, Tessier et al. 1993, Lovejoy et al. 2001). The UMM is a generalization of the lognormal model (e.g. Kolmogorov 1962, Monin and Yaglom 1975) and relies on using the Levy family distribution. The Levy distribution depends on four parameters, and one may write a Levy noise as S (α, C1, µ, β) (Janicki and Weron 1994). For multifractal applications, the centering term µ=0 and the skewness parameter β=-1 (Schertzer and Lovejoy 1989). The value adopted for this latter parameter results in a distribution extremely skewed to the left. Therefore, a generated Levy noise for multifractal applications depends only on two parameters: the Lévy index α and the intermittency parameter C1, which can be estimated from the following equation: 1.11 𝜉 𝑞 = 𝑞𝐻 − !!!!! 𝑞! − 𝑞 , 𝛼 ≠ 1 where ζ(q) is the scaling exponent defined for any statistical moment q. C1, usually known as the scale parameter (0<C1<1), is equal to half the variance when α=2 and also plays a similar function when α<2. In other words, it measures the width of the distribution. As C1 increases the magnitude of the sudden large jumps increase (Boufadel et al. 2000), hence a value of C1 near 1 corresponds to a process that usually has very small observations, except 12 on rare and distant occasions when an observation significantly exceeds the average (Seuront and Lagadeuc 2001). The parameter α, also known as the Lévy index, measures the degree of multifractality and determines the shape of the probability distribution (Schmitt et al. 1998). A value of α=0 correspond to the monofractal case and α=2 (its upper limit) to the lognormal multifractal case (Tessier et al. 1993). As α decreases, the frequency of sudden large jumps in the random field increases (Tennekoon et al. 2003). In particular, for the monofractal case, that is, when C1=0 or α=0, ζ(q) is a linear function of q. 1.5 Fractal analysis in marine ecology There are several studies which have focused on fractal analysis in marine ecological research. A good text, which summarises many of them, is provided by Seuront (2010). Here we briefly describe some of them which focus on benthic and planktonic populations as examples of fractal research studies conducted during the last three decades. Bradbury et al. (1984) analyzed the hierarchical organization of a coral reef in Australia using the fractal dimension D to represent transition zones in the spatial organization of the reef. Azovsky et al. (2000) studied the spatial distribution of intertidal communities of macrozoobenthos and microphytobenthos in Kandalaksha Bay (White Sea, Russia). They found that these communities have a fractal spatial structure that can be described as a hierarchy of nested patches of various sizes from meters to kilometers. In a similar study, Seuront and Spilmont (2002) analyzed the spatial distribution of the biomass of several genera of microphytobenthic diatoms from the intertidal zone of an exposed sandy beach in Wimereux, France (Eastern English Channel). They found that the distribution of the biomass of these diatoms showed a very intermittent behavior, which can be characterized 13 by a few dense patches and a wide range of low-density patches of different size, where the frequency of occurrence of patches of different concentrations is related to the size of patches by a fractal dimension D. The existence of this self-similar behavior is, according to Seuront and Spilmont (2002), indicative of self-organization near a phase transition, an ecological state where large-scale correlations (i.e., long-range dependence) can emerge. Seuront and Lagadeuc (2001) analyzed the variability in the temporal abundance of the calanoid copepod Temora longicornis sampled from the coastal waters of the Eastern English Channel (France). They estimated the fractal dimension of the abundance of T. longicornis using the power spectrum analysis (D=1.79), and estimated the universal multifractal parameters (H, α, C1) using the structure function approach. They found a single scaling regime for a range of statistical moments q from 0 to 3, which indicates that the same process can be regarded as being responsible for the scaling of the abundance of this species from 6 minutes to 66 hours. The reported scaling exponent for the first statistical moment ζ(1), was equal to 0.58, which indicates that the temporal distribution of T. longicornis was not stationary. They also found that the fractal dimension of T. longicornis (D=1.79) is: (i) significantly higher than those estimated for oceanic turbulence (D varies from 1.2 to 1.4; Osborne et al. 1989, Osborne and Caponio 1990), (ii) higher than those of phytoplankton distributions estimated from in vivo fluorescence time series in the Southern Bight of the North Sea and the Eastern English Channel (1.61<D<1.67; Seuront et al. 1996a,b; Seuront et al. 1999), and (iii) very similar to that estimated for the abundance on transects of the copepod Neocalanus cristatus from the subarctic Pacific over a similar range of scales (D=1.80, Tsuda 1995). They interpreted these results by suggesting that differences in size and motility between phytoplankton and zooplankton can affect copepod 14 behavior (e.g. diel migration, phototaxis, social behavior, predation pressure, etc.) inducing larger fractal dimension (i.e., a large degree of occupation of a tri-dimensional Euclidian space) for zooplankton species. The similar D’s obtained for T. longicornis and N. cristatus suggest that the distribution of both zooplankton species could be independent of their surrounding environment. On the other hand, the universal multifractal parameters estimated in this study (H=0.58 and C1= 0.44) were higher than those of phytoplankton biomass (0.12<H<0.34, 0.02< C1<0.05, Seuront et al. 1996a,b; Seuront et al. 1999) implying that zooplankton distributions are less stationary and sparser. The high degree of α (1.60) indicates that intermittencies of all magnitudes contribute to the multifractal abundance dynamics of T. longicornis. Seuront and Lagadeuc (2001) recognized the usefulness of the universal multifractal approach, which does not require any preconceptions of the data. In this way, it can be considered as a real alternative to quantify intermittency and non-Gaussian statistics derived from patchy distributions frequently found in plankton communities. 1.6 Spatial fractal analysis in fisheries oceanography Spatial multifractal dynamics for a predator-prey relationship that varies as a function of the predation rate has been reported in simulated fish populations (Tikhonov et al. 2001). In oceanic ecosystems, fishing vessels performed Levy flights, a type of motion characterized by long-tailed distributions with no characteristic scale for the observed distribution, i.e. it is scale invariant (Bertrand et al. 2005). Fish shoals of different fish species follow a power-law or fractal spatial distribution, which means that an instantaneous structural similarity exists at all scales observed, from tens of meters to tens of kilometers. Similar 15 underlying behavioural mechanisms at all scales could be the cause of this spatial configuration. In other words, this scale invariance relationship means that very large shoals are structurally similar to much smaller fish groups (Makris et al. 2006). McClatchie et al. (1994) investigated the spatial variability in biomass of Antarctic krill (Euphausia superba) estimated in four acoustic surveys conducted around Elephant Island in the Southern Ocean. Using fractal geometry, these researchers estimated the Haussdorf dimension9 of the spatial variability of abundance estimates in order to obtain new information about the appropriate sampling scales to analyze the spatial patchiness of krill swarms. In a similar study, Wang et al. (1999) reported a fractal spatial variability of the catch rates obtained from the Northern Prawn Fishery, Gulf of Carpentaria, Australia. These authors hypothesized that higher catch rates were related to the complexity of the reef (untrawlable grounds) rather than the total area of the fishery. To validate their hypothesis, they estimated fractal dimensions from different areas (as a measure of complexity of those areas) to explain the variation of catches among areas. Preliminary results indicated that differences between fractal dimensions could be associated with the complexity of the fishing grounds and also with the catch rates. Snover and Commito (1998) studied the spatial complexity and distribution of an extensive blue mussel (Mytilus edulis) bed in a population from an intertidal flat of sandy mud in Maine (USA) taking areal-view photographs above the bed. Using quadrats that were placed in the field, they estimated the fractal dimension (D) of the mussel spatial patterns. The spatial pattern of mussels in every quadrat was fractal with scaling ranges over nearly two orders of magnitude. Furthermore, the fractal dimensions were significantly related to 9 The Haussdorf dimension is a measure of how similarly distributed a variable is when viewed at different levels of resolution. 16 the mussel percentage cover and densities within the quadrats. These relationships were curvilinear (i.e., dome shaped) with D ranging from 1.36 to 1.86. The highest fractal dimensions were found at intermediate levels of both percentage cover (approximately 50%) and densities (approx. 240-250 ind/0.0625 m2). They also found the fractal dimension (D) decreased as spatial aggregation increased. Snover and Commito (1998) suggested that the fragmented patch structure, which relates to a higher fractal dimension of the mussel bed, might be caused by physical disturbances such as storms and ice scouring. Guichard et al. (2003) analyzed the spatial dynamics of another species of mussel (Mytilus californianus) from the Oregon coast using theoretical and experimental data to describe and quantify the disturbance and recovery of the mussel beds. Mussels are disturbed by strong waves, especially during winter, that create gaps in the mussel bed. The model that simulates the theoretical data included local interactions among elements of the mussel bed (i.e., biotic processes) and between the mussel bed and wave disturbance (i.e., abiotic processes). They found, for simulated and field data sets, clear scale invariant patterns in the mussel beds characterized by power law distributions. Τhese fractal patterns, characteristic of systems with long-range correlations, emerge as a response to oceanographic conditions that allow self-organization in benthic communities. 1.7 Temporal fractal analysis in fisheries oceanography As discussed in previous sections, fractal analysis has been focused on the study of the spatial structure of marine populations including fish resources of economic interest. On the other hand, the application of fractals to the study of marine population dynamics has been strongly focused on planktonic time series, and many of them conducted by Dr. Laurent 17 Seuront and colleagues (see Seuront 2010 and references therein). The study of fisheries time series using a fractal approach has evolved slowly and, as only a few applications were found in the scientific literature, we describe them with more detail in the following paragraphs. 1.7.1 Case study 1 (Ferriere and Cazelles 1999) The alternation of periods of resource commonness and periods of scarcity is a well-documented phenomenon in several commercial fisheries (Soutar and Isaacs 1974, DeVries and Pearcy 1982, Rothschild 1986). The intermittent dynamics of the abundance patterns of some fish captured in fisheries are qualitatively similar when observed at different time scales (e.g., the abundance pattern of Sardinops caerulea over 12,000 years is similar to the pattern observed over the past 2000 years) which give rise to the analysis of the temporal dynamics of intermittent rarity10 of fish populations. In this context, Ferriere and Cazelles (1999) using simulated and real data of the Pacific sardine (Sardinops caerulea) found a well-defined scaling relationship (exponent equal -3/2) in the distribution of time spent in phases of rarity. This power law implies that there are rarity phases of arbitrary length (i.e., the statistical distribution of rarity times is scale-invariant), and also that the temporal patterns are similar at all density levels, which suggest the presence of a fractal pattern in intermittent rarity time series. 1.7.2 Case study 2 (Montes 2004) The temporal dynamics of catch (C) and standardized catch per unit effort (CPUE) of the 10 Intermittent rarity corresponds to the alternation over variable periods of time phases of extremely low abundance and short outbreaks. In this context, rarity times correspond to the sequence of times elapsed between successive outbreak events. 18 southern Hake (Merluccius australis) fishery off Chile was analyzed by Montes (2004) from a fractal perspective. Taking into account the dynamics of daily catch rates and effort, and also the reduced fishing effort due to the passage of the Chilean General Fishing and Aquaculture Law in 1991, the Hurst coefficient (H) and fractal dimension (D) of daily C and CPUE time series from different stages (periods) in the history of the fishery were estimated. Using Dispersional Analysis (Caccia et al. 1997, Cannon et al. 1997) they found significant changes in the magnitude of these parameters between time periods (e.g. catch time series: D=1.20, H=0.80 for the 1980-1984 period; D=1.34, H=0.66 for the 1985-1989 period; D=1.28, H=0.72 for the 1990-1994 period). On the other hand, using structure functions and the Universal Multifractal Model (UMM, Schertzer and Lovejoy 1987, equation 1.11) multiple well-defined scaling relationships between 2 and 256 days (2.11 decades, equation 1.3) for the daily C and between 2 and 64 days for the daily CPUE time series for several statistical moments (i.e., from 0 to 7) were detected. The UMM used by Montes (2004) to quantify the dynamics of the southern hake time series fully incorporates intermediate scales of temporal variability (e.g. from weeks to seasons) which can play a fundamental role in the dynamics of this exploited population. The main hypothesis behind the results of this study was that the observed multifractal dynamics for the catch time series reflect multifractal predator-prey temporal dynamics between the trawl factory fleet and the southern hake that mainly depend on the predation rate of the fleet. In this context, the use of fractal and multifractal parameters as indicators of change in the status of this fishery was proposed. For example, the relationship between the scaling exponent ζ(q), the degree of multifractality (α), and the standard deviation (SD) of southern hake C time series to discriminate between catches of time periods with 19 different fleet-resource (predator-prey) dynamics was used. 1.7.3 Case study 3 (Halley and Stergiou 2005) Halley and Stergiou (2005) measured the variance of fisheries landings using 344 time series extracted from the Food and Agriculture Organization of the United Nations (FAO) database for the period 1950-1996 and compared their results (i.e., the observed variance growth in landings) with those obtained using 431 terrestrial population time series from the Global Population Dynamics Database (GPDD; Inchausti and Halley 2001, 2002). In order to assess the rate of increase of the variance as a function of the length of the time series they estimated the Hurst coefficient (H) of marine and terrestrial populations. For long-range dependent (LRD) time series, the variance (V) grows with the time interval Δt (Feder 1988) as: 1.12 𝑉 ∝ |∆𝑡|!! where the exponent 2H (and also H) could be estimated from the following equation: 1.13 log𝑉! = 2𝐻 𝑙𝑜𝑔 𝑘 + 𝑐 where 𝑉! is the average of time series variance calculated using subsequences of size k (time series resolution). For the majority of the stocks analyzed (e.g., Limanda limanda, common dab in NE Atlantic; Trachurus murphyi, Chilean jack mackerel in SE Pacific; Lophius americanus, American angler in NW Atlantic) landings variability increased with the length of time over which it was calculated. For example, the number of stocks characterized by a variance of 0.5 or higher increased from 6 (estimated for a three year period) to 64 20 (estimated for a 47 year period). This significant proportion of H values greater than 0.5 (implying an accelerating variance growth) is associated with a strong reddening of the landings data. However, these results must be analyzed very carefully because for some landings time series (i.e., jack mackerel), after removing the effect of an exponentially growing fishery (which causes a linear increase in landings on a logarithmic scale associated with an increasing fishing mortality), the H coefficient went below 0.5 (implying a decelerating variance growth). On the other hand, Halley and Stergiou (2005) conducted a similar analysis using the scale deposition rates of Pacific sardine and Northern anchovy (data obtained from Baumgartner et al. 1992), which showed a steady rise in variance as a function of the observation time to scales of centuries and above. Halley and Stergiou (2005) recognized that the increasing variability in the landings data represents a significant problem for fisheries management because it will introduce an additional level of uncertainty into various estimations. Indeed, most traditional fisheries management models presently used assume the existence of an equilibrium yield, which could not be reached according to the unbounded variance growth of the landings time series estimated in this research. 1.7.4 Case study 4 (Niwa 2006, 2007) Long-range correlations in annual recruitment time series were analyzed by Niwa (2006) using Rescaled-Range Analysis. Hurst exponent was calculated for the ICES-ACFM (International Council for the Exploration of the Sea-Advisory Committee on Fishery Management) time series data that covers 27 fish stocks from the North Atlantic Ocean, including: (i) Cod (Gadus morhua), (ii) Greenland Halibut (Reinhardtius hippoglossoides), (iii) Haddock 21 (Melanogrammus aeglefinus), (iv) Herring (Clupea harengus), (v) Plaice (Pleuronectes platessa), (vi) Saithe (Pollachius virens), (vii) Sole (Solea vulgaris), (viii) Sprat (Sprattus sprattus) and (ix) Whiting (Merlangius merlangus). This author found a Hurst coefficient equal to 0.87 for all averaged stocks, which implies the existence of long-range correlations for fish recruitment time series. The same author analyzed the variability in population growth rate from the same group of North Atlantic fish stocks but did not find any evidence of correlation among years (Niwa 2007). He concluded that changes in population trajectories are independent from one year to the next, and can be described by a geometric random walk process. 1.7.5 Case study 5 (Mendes et al. 2014) The same database from the North Atlantic Ocean analyzed by Niwa (2006, 2007) was used by Mendes et al. (2014) to re-examine the dynamics of population growth rates of three fish stocks. In this case, fish time series for each stock were analyzed separately, in contrast to the averaging procedure of time series conducted by Niwa (2006, 2007). Hurst coefficients significantly greater than 0.5 were found for Herring, Cod and Haddock stocks (H fluctuated between 0.8-0.9) which means that the geometric random walk hypothesis (i.e. the growth rate is independent from one year to the next) proposed by Niwa (2006) was rejected. The long-range dependence (LRD) dynamics for these three stocks could be associated with different exploitation regimes and also with large-scale changes in environmental conditions produced by the North Atlantic Oscillation which can induce low or high productivity regimes in fish recruitment (Mendes et al. 2014). 22 Under the framework of the fractal theory summarized in previous paragraphs, and according to the evidence of fractal temporal patterns of variability found in fishery time series worldwide, the following four main objectives were proposed. 1.8 General objectives 1. To determine and quantify fractal patterns of variability of the smooth pink shrimp fishery off the west coast of Vancouver Island. 2. To determine and quantify fractal patterns of variability of oceanographic time series. 3. To determine and quantify bivariate fractal or multifractal relationships between oceanographic and smooth pink shrimp fisheries time series. 4. To evaluate the usefulness of fractal parameters as indicators of change in the dynamics of Pandalus jordani fishery. For testing Objective 1, a three fishing season time series of daily catches of the smooth pink shrimp (Pandalus jordani) fishery was used in Chapter 2. Then in Chapter 3, a more extended time series of catches was used to confirm previous results and to detect bivariate fractal or multifractal relationships with sea surface temperature and wind stress time series from the west coast off Vancouver Island (Objectives 2 and 3). Finally, in Chapter 4 a fractal-based indicator was constructed and tested against other indicators. To accomplish Objective 4, the temporal dynamics of this latter indicator was analyzed in relation to changes in the Pandalus jordani fishery and proposed to detect early changes for this fishery. Chapter 5 summarizes previous results and gives recommendations for future studies in fractal research for this and other fisheries. 23 Chapter 2: Multifractal patterns in the daily catch time series of smooth pink shrimp (Pandalus jordani) from the west coast of Vancouver Island, Canada11 2.1 Introduction Fractal theory has been used in aquatic sciences to quantify scale-invariant relationships under the form of scaling or power laws (Seuront 2010). These scaling laws represent patterns that help to understand the complex structure and organization of marine ecosystems and their populations through space or time (Pascual et al. 1995, Schmid 2000). Under this framework, scale-invariant patterns are detected in a time series by quantifying the dependence between observations and integrating many different temporal scales into a single model. This approach has been proven to be useful in the identification and modelling of ecological and environmental processes that affect the temporal variability of marine planktonic populations (Lovejoy et al. 2001, Fisher et al. 2004). In fisheries research, the most common approximations in the statistical modelling of time series involve: (i) modelling of short-range correlations using univariate (autoregressive integrated moving average, ARIMA) and multivariate (transfer function noise, TFN) models (e.g. Downton and Miller 1998, Hanson et al. 2006), and (ii) estimation of short-range correlations (or calculation of the number of independent observations) (e.g. Pyper and Peterman 1998; Perry et al. 2000). The first approach can detect memory patterns at 11 A version of this chapter was published in Canadian Journal of Fisheries and Aquatic Sciences as: Montes, R.M, Perry, R.I., Pakhomov, E.A., Edwards, A.M. and J.A. Boutillier. 2012. Multifractal patterns in the daily catch time series of smooth pink shrimp (Pandalus jordani) from the west coast of Vancouver Island, Canada. Canadian Journal of Fisheries and Aquatic Sciences 69: 398-413. 24 short time scales and can use those patterns to forecast future values, whereas the second approach is often used to adjust the degrees of freedom after conducting classical linear regression and correlation analyses. Both approximations are designed to detect the correlation or dependence between observations that are closely located to each other (i.e., short-range correlations), ignoring the possible dependence that could exist between distant observations (i.e., long-range correlations). Long-range correlations have been detected in fisheries time series (e.g. Halley and Stergiou 2005) but the high variability (intermittency) and non-linearity usually found in marine catch time series have not been taken into account. These complex dynamics can be quantified in fisheries time series using a framework of long-range correlations developed for studies of turbulence and applied to plankton research (e.g. Seuront et al. 1999, Seuront and Lagadeuc 2001). Many attempts have been made to detect and model long-range correlations in the field of marine ecology, such as in characterizing the spatial distribution of exploited and non-exploited marine populations and communities (e.g. Azovsky et al. 2000; Guichard et al. 2003). Fractals have been used extensively to model long-range correlations and the intermittent dynamics of marine planktonic and oceanographic time series (e.g. Seuront et al. 1999; Lovejoy et al. 2001; Fisher et al. 2004). However, despite the ability of fractals to integrate multiple time scales of observation in a given model, few examples of their use in modelling marine exploited populations in the time domain can be found in the fisheries literature. Two of these examples are the studies by Halley and Stergiou (2005) and Niwa (2007) who focused on analyzing the dynamics of annual patterns in fishery landings and marine population growth rates, respectively. In a similar way, Niwa (2006) reported long-range correlations in annual recruitment time series of fish populations from the North- 25 Atlantic. Since these studies were published, better methods for estimating the degree of long-range correlations for time series on the order of 100 data points (Chamoli et al. 2007) and shorter (e.g. 64 data points, Delignieres et al. 2006) have become available. In addition, methods are now available to estimate how far the effect of long-range correlations are detectable (i.e. extension of scaling ranges) (e.g. Seuront et al. 2004) quantitatively rather than visually as has been done in previous examples from the fisheries literature. In the present study, we analyzed the dynamics of fisheries catch time series using a finer scale of temporal resolution (daily) than previous fisheries studies, testing for the existence of multiple scaling patterns (mutifractality). We use data on daily catches from the commercial smooth pink shrimp (Pandalus jordani) fishery off the west coast of Vancouver Island, Canada, for the period 1994-1996. These data can be considered as long (N=642) because the daily resolution allows for the analysis of intra-annual scales of variability. As the main objective of this thesis is to construct an early warning indicator based on fractal geometry, first the fractality of smooth pink shrimp fishery catches should be established. In this context, the first goal of this Chapter was to demonstrate the existence of fractal properties in daily fisheries catch time series using the longest continuous segment of daily catches that extends from 1994 to 1996. Our second goal focused on comparing the monofractal model with the more complex multifractal model regarding their abilities to detect and quantify long-range correlations and high variability. A multifractal model was incorporated into our analysis because it allows the detection of non-linear and intermittent fluctuations usually found in fisheries time series, in addition to the detection of long-range correlations. Our data come from the commercial smooth pink 26 shrimp (Pandalus jordani) fishery off the west coast of Vancouver Island, British Columbia. This fishery was generally open year-round up to 1997 with no quotas from the late 1980s to 1997. The fishery was conducted by small vessels using otter and beam trawls, with footrope lengths between 9 and 27 m and vertical openings greater than 1.5 m (Perry et al. 2000; Rutherford et al. 2004). Pandalus jordani is a protandric hermaphrodite species, which generally lives up to 3 years; rarely up to 4 years. This means that they spawn first as males and then change their sex to females around 2 years of age (Butler 1980). Spawning occurs during the fall-winter period and the eggs hatch during the spring (Butler 1980). Individuals fully recruit to the fishery at the age of 2 years (Boutillier et al. 1997). We start with a basic introduction to fractals and long-memory processes, and then analyze the pink shrimp trawl fishery data for the presence of these long-range correlations and intermittent dynamics. We conclude by discussing some hypotheses that may explain the observed patterns of variability. We also discuss the potential for this approach to be used to assist with the management of this fishery and to develop early warning indicators of changes in the underlying environmental or human processes within which this fishery operates. 2.2 Methods 2.2.1 Fractals and long-memory processes A fractal can be defined as an object that lacks a characteristic length because the same geometrical structures are present at all length scales. If the observational scale is reduced, new structures arise, and the new ones are similar to those of the original unit. As a consequence, fractals are composed of subunits (and further subunits) that resemble the 27 larger scale structure, a property known as self-similarity (Mandelbrot 1983, Goldberger et al. 2002). Self-similar processes are invariant in distribution under scaling of space and/or time, and the scaling coefficient (or index of self-similarity) is a positive number denoted H (the Hurst coefficient; Samorodnitsky and Taqqu 1994). Here, 0<H≤1. Due to this property, the shape of fractals is nonrectifiable, meaning that they consist of an infinite sequence of clusters within clusters or waves within waves. In rectifiable objects, increasingly accurate measurements based upon successively smaller scales converge to a limit which is the true extent (length) of the object (Mandelbrot 1983). However, in fractals, the same procedure generates an infinite series according to a power-law: 2.1 𝜆 𝑟 = 𝑘𝑟!!! where λ is the length of the object measured at the scale unit r (the length of the object diverges as r → 0), k is a constant, and D is the fractal dimension. The definitions of all variables are summarized in Table 2.1. D exceeds the topological dimension (d) of an object and is not an integer, and satisfies d<D<d+1, where d+1 is the space dimension of the object. For example, for a fractal time series 1 <D< 2, and D=2-H. 28 Table 2.1. Definitions of variables and parameters used in Chapter 2. Symbol Definitions H Hurst coefficient C1 Intermittency parameter α Lévy index λ Length of a fractal object r Time lag or scale unit (fishing days) k Constant of power law equation D Fractal dimension d Topological dimension (d =1 for time series) ui Total daily catch observations of the linearly detrended time series n, m Sample size and number of cosine curves, respectively a0 ; aj, bj Intercept and amplitudes of Fourier series cj Periods of Fourier series i,j Daily values and Fourier index, respectively w Frequency 29 Symbol Definitions q Statistical moment β Power spectrum slope xi Value of the TDCa time series for day i Yr(q) Structure function empirical observations Ŷr(q) Structure function fitted values N Total length of pink shrimp catch time series ζ(q) Structure function scaling exponent H(q) Hierarchy of Hurst coefficients Κ(q) Scaling moment function U Theil’s inequality coefficient r1 Lower scaling limit r2 Upper scaling limit R’ Deviance p1, p2 Number of parameters of monofractal and multifractal models α’ Significance level 30 Fractal (or scale-invariant) processes have no characteristic scale, which means that all scales contribute to the observed dynamics. This property is referred as a scaling law or scaling behavior (Abry et al. 2009). Fractal processes generate irregular fluctuations on multiple time scales, analogous to fractal objects that have a wrinkle structure on different length scales (Mandelbrot 1983). In an idealized model, scale-invariance holds on all scales, but the real world imposes upper and lower bounds over which such behavior applies (Goldberger et al. 2002). In contrast to monofractals, which are homogeneous in the sense that they have the same scaling properties characterized by only one exponent or coefficient (i.e., H) throughout the entire signal, multifractals are a class of signals which requires a large number of indices to characterize their scaling properties. Multifractal signals are intrinsically more complex and inhomogeneous than monofractals, i.e., different parts of the signal may have different scaling properties (Goldberger et al. 2002). Monofractal analysis looks at the geometry of a pattern, whereas multifractal analysis looks at the arrangements of quantities. A monofractal can be defined as a geometrical set of points and a multifractal is defined as a mathematical measure (Lavallée et al. 1993). The multifractal approach implies that a statistically self-similar measure can be represented as a combination of interwoven monofractal sets with corresponding scaling exponents. A combination of all the monofractal sets produces a multifractal spectrum that characterizes variability and heterogeneity of the analyzed variable (Kravchenko et al. 1999). In this study, to determine whether a monofractal or a multifractal model better describes pink shrimp daily catches, we used an objective statistical procedure that has not been applied previously in this context. 31 One means of capturing long-range correlations or long-memory features of fractal time series is with self-similar stochastic processes, also known as stochastic fractals (Mandelbrot 1983, Samorodnitsky and Taqqu 1994). A long memory process has an autocovariance function with a decay rate much slower than the decay rate characteristic of a standard stationary autoregressive process (Percival et al. 2001). In other words, the observations separated by long time intervals still exhibit a non-zero covariance (current observations retain some “memory” of distant past observations). Self-similar models have been studied for more than three decades and have been successfully applied in a variety of fields in the form of fractional Gaussian noise (fGn) and fractional Brownian motion (fBm) (e.g., Mandelbrot and Van Ness 1968, Rodriguez-Iturbe and Rinaldo 2001). For self-similar time series the Hurst coefficient (H) is a measure of long-range dependence. When H>0.5, all the observations of a time series are positively correlated, and the closer H is to 1 the smoother the function. When H<0.5, all observations are negatively correlated. A time series having the characteristics of Gaussian white noise corresponds to a fractional Gaussian noise series with H=0.5 (Samorodnitsky and Taqqu 1994). Fractional Gaussian noise corresponds to a stationary process with constant mean and variance, whereas fractional Brownian motion is non-stationary with stationary increments (Mandelbrot and Van Ness 1968). Differencing fBm creates fGn, and cumulatively summing fGn produces fBm. Both processes are characterized by the same Hurst coefficient but display different dynamics. 2.2.2 Pink shrimp fishery and total daily catches The smooth pink shrimp (Pandalus jordani) is captured by the trawl fishery with six other species of shrimp and, up to 1996, the majority of the catch was a mix of smooth pink 32 shrimp (greater than 90%) and sidestripe shrimp (Pandalopsis dispar) (Rutherford et al. 2004). Individual trawls (weight of shrimp in kilograms for each trawl) for P. jordani were obtained from the Department of Fisheries and Oceans, Canada. This time series takes into account Otter and Beam trawls and is from offshore Pacific Fisheries Management Areas (PFMA) 121, 123, 124 and 125 (Fig. 2.1). Catch data were recorded by fishermen in harvest logbooks at the time of fishing. Figure 2.1. Fisheries and Oceans Canada Pacific Fishery Management Areas (PFMA). Data for this study came from PFMA 121, 123, 124 and 125 off the west coast of Vancouver Island, Canada. 33 A total daily catch (TDC: the daily sum of recorded catches) time series was constructed for the west coast Vancouver Island smooth pink shrimp fishery for the 1994-1996 period. Only individual trawls registered during the principal spring/summer fishing season (April-October, Perry et al. 2000) were used, giving 214 days each year. No major fishing activity was recorded during the November-March period of each year. In consequence, sporadic days with trawling activity during this latter period were not considered as fishing days and were not included in this study. Total daily catches for the three years were ‘stitched’ together, to yield a total of 642 fishing days (Fig. 2.2); the consequences of the stitching will be discussed later. Figure 2.2 Pink shrimp (Pandalus jordani) total daily catch (i.e. daily sum of all recorded catches; TDC) time series from 1994 to 1996 for the April-October period (642 fishing days) in Pacific Fishery Management Areas 121, 123, 124, 125 west of Vancouver Island, Canada. A sinusoidal curve with an annual period of 214 fishing days was fitted to this time series. A slightly increasing linear trend is also visible in this time series. 34 This time series corresponds to the longest period of total daily catches from the main fishing season for which the smallest number of days with no observations occurred (24 “missing” observations, equivalent to 3.7% of the total). Missing observations were estimated by linear interpolation using one point on either side of the missing value. Interpolations can add correlations to time lags less than 30 data points for time series with 10% of missing data filled by linear interpolations (Wilson et al. 2003). In our case, as the fraction of missing observations is less than 5%, we expect to find the effect of spurious correlations (if any are present) at scales significantly less than one month. The pink shrimp TDC time series presented a marked periodicity and a slight increasing linear trend (Fig. 2.2). Its periodogram (not shown) has a strong peak at a frequency (0.0047 fishing days-1) equivalent to the period of one fishing season (214 fishing days). This corresponds to an annual cycle (summer to summer) due to the stitching together of the three years of data. Periodicities are considered as deterministic trends, which can obscure the detection of scaling behaviour, and as such they must be removed before the estimation of long-memory parameters (Montanari et al. 1999). First, the linear trend of the TDC time series was removed by a least squares fit (R2=4%, p<0.01). The seasonal period corresponding to the annual cycle (frequency equal to 0.0047 fishing days-1) and its first harmonic (frequency equal to 0.0093 fishing days-1), which in total explained 25% of the time series variance, were identified in the resulting TDC linearly detrended time series by fitting a Fourier series with m cosine curves using the following equation: 2.2 𝑢! = 𝑎! + 𝑎!𝑐𝑜𝑠 2𝜋𝑖𝑗𝑛 + 𝑏!𝑠𝑖𝑛 2𝜋𝑖𝑗𝑛!!!! ; 𝑎! = 𝑏! = 2𝑁 𝑢!𝑐𝑜𝑠 2𝜋𝑖𝑗𝑁!!!! 35 Here, cos(2πij/n) and sin(2πij/n) are used, for a fixed frequency, as predictor variables to fit the intercept a0 and the amplitudes aj and bj from the data ui using an ordinary least squares regression method (Cryer and Chan 2008). The time series of residuals obtained after the elimination of deterministic trends is hereafter called the TDC residual time series (TDCr). A step-by-step diagram showing all transformations conducted on the pink shrimp TDC time series is shown in Fig. 2.3. The above linear detrending procedure and the removal of seasonal components are Step 1. For Step 2, nonstationary segments (time intervals showing a marked trend) in the TDCr time series were removed as they tend to give a biased estimate of H near 1 (Yu et al. 2003). After that, the mean of TDCr was subtracted from the individual values and the cumulative sum (partial summation) of the obtained values was calculated (e.g. Gao et al. 2006). The series obtained after this operation is hereafter called the TDC anomaly time series (TDCa). The TDCa time series is free of periodicities and is non-stationary, as is required to estimate the multifractal parameters. 36 Figure 2.3. Diagram showing the transformations and mathematical operations conducted on pink shrimp (Pandalus jordani) total daily catch time series (TDC) on a step by step basis. TDCr represents the total daily catch residual time series and TDCa is the total daily catch anomaly time series obtained after cumulatively summing the TDCr time series. 37 2.2.3 Monofractal and multifractal analysis Scaling dynamics and stationarity were tested using power spectrum analysis (Davis et al. 1994, Rodriguez-Iturbe and Rinaldo 2001). If scaling is present in pink shrimp catches, the characteristic power-law form of fractals will appear as: 2.3 𝐸 𝑤 ≈ 𝑤!! where w is frequency, E(w) is the power of the TDCa time series (squared modulus of the Fast Fourier Transform; Davis et al. 1996) and β an exponent to be estimated. On a logarithmic scale, scaling dynamics will manifest as an approximately linear curve when plotting E(w) versus w. Power spectrum analysis is equivalent to an analysis of variance in which the total variance (statistical moment q=2) is partitioned into contributions coming from processes with different time scales giving a quantitative measure of the importance of each frequency (Seuront et al. 1999). As a consequence, highly intermittent and non-Gaussian time series are not appropriately quantified using this methodology, which focuses only on the second order moment (Seuront and Lagadeuc 2001). To completely quantify the scaling dynamics of the TDCa time series we analyzed higher statistical moments using a function of order q (statistical moment) known as a structure function (Davis et al. 1996, Tennekoon et al. 2005) or generalized variogram (Lavallée et al. 1993, Rodriguez-Iturbe and Rinaldo 2001). When q=2 the structure function corresponds to the variogram (Davis et al. 1994). Structure functions are only applicable to non-stationary self-similar time series (Davis et al. 1994, Yu et al. 2003), i.e., a time series for which the slope β obtained from equation (2.3) fluctuates between 1<β<3. These time series correspond to nonstationary processes with stationary increments (Davis et al. 1994, 38 1996). However, if -1<β<1 the time series is classified as strictly stationary and a transformation to convert it to a nonstationary process is necessary before using structure functions (Yu et al. 2003, Gao et al. 2006). These functions have been used in plankton ecology (e.g. Seuront et al. 1999, Fisher et al. 2004) and were also applied in the analysis of the multifractal temporal dynamics of the southern hake fishery of Chile (Montes 2004). The structure function Yr(q), for time lag r and statistical moment q, is defined as: 2.4 𝑌!(𝑞) = 𝑥!!! − (𝑥!) ! ; 𝑟 = 1, 2, 3,… ,𝑁 4 where N is the length of the TDCa time series, xi is the value of the TDCa time series for day i, and the operator 〈...〉 corresponds to the mean over i = 1,2,3,…,N-r. The time lag r is equivalent to the scale unit defined in equation (2.1) and is measured using previously defined fishing days. The maximum time lag r is usually set to N/2 because a greater time lag is more affected by low sample sizes and also by spurious properties of the data (Journel and Huijbregts 1978). In order to be conservative, in this study the maximum time lag r was set to N/4. In the averaging procedure, we considered 482 pairs of data from the TDC anomaly time series, which is significantly more than the 100 data pairs (Tennekoon et al. 2005) or 225 data pairs (Webster and Oliver 1992) considered reliable when structure functions or their equivalent variograms are calculated. The scaling exponents ζ(q) were estimated by the slope of the linear regression between Yr(q) and r on a logarithmic scale for different statistical moments q using the following equation: 2.5 𝑌!(𝑞) ≈ 𝑟! ! 39 We considered 40 moments in the structure function calculations, such that q varies between 0.1 (weakest and most frequent fluctuations) and 4.0 (strongest but sporadic fluctuations) with an increment of 0.1. We then define a hierarchy of Hurst coefficients using the ζ(q) values as: 2.6 𝐻 𝑞 = ζ(𝑞)𝑞 which is the goal of this method (Davis et al. 1994). Those processes with constant Η(q), characterized by simple scaling, are known as “monofractals” (or simple fractals). Processes with nonlinear and convex (downward facing) ζ(q) are called “multifractals” (Seuront et al. 1999, Tennekoon et al. 2005). For monofractals, ζ(q) is linear (e.g. ζ(q) = q/2 for Brownian motion). In particular, the first moment (q=1) gives the scaling exponent H = ζ(1), corresponding to the scale dependency of the average fluctuations. If H ≠ 0 the fluctuations Yr(q) will depend on the time scale, which characterizes the degree of stationarity of the process (Seuront and Lagadeuc 2001). No scaling is present for a strictly stationary process ζ(q)~0 and therefore structure functions are also used to detect stationary (and nonstationary) segments in a time series before the estimation of multifractal parameters (Yu et al. 2003). We used the Universal Multifractal Model (Schertzer and Lovejoy 1987, Lovejoy et al. 2001) to characterize and quantify the dynamics of a multifractal process. In using this model we avoided estimating a large number of exponents, ζ(q), to completely characterize the dynamics of the pink shrimp catch time series. Under this framework, the multiscaling 40 behaviour of an intermittent process can be quantified using a scaling moment function K(q) which depends only on two parameters, C1 and α, estimated through: 2.7 𝐾 𝑞 = 𝐶!𝛼 − 1 𝑞! − 𝑞 ; 𝛼 ≠ 1 𝐶! 𝑞 log 𝑞 ; 𝛼 = 1 The intermittency parameter C1 characterizes the heterogeneity of the process; when this parameter increases, the magnitude of sudden large jumps increases (Seuront and Lagadeuc 2001). It is also known as the codimension, i.e. the fractal dimension is d - C1 (d=1 for time series). A value of C1 near 0 indicates a process (in this case, a time series) having near-average observations at all times (homogenous), whereas a value of C1 near 1 corresponds to a process that usually has very small observations except on rare and distant occasions when an observation significantly exceeds the average (i.e., heterogeneous or sparse) (Seuront and Lagadeuc 2001). It goes beyond the standard concept of asymmetry known as skewness (which only looks at the third statistical moment) and has a long-term memory signature (Seuront 2010). The parameter α, also known as the Lévy index, measures the degree of multifractality and determines the shape of the probability distribution; it satisfies 0<α≤2 (Seuront et al. 1999). As α decreases, the frequency of sudden large jumps increases. For a multifractal, the departure from the simple linear behaviour can be quantified by extending equation (2.7) to include the scaling function K(q) as: 2.8 ζ 𝑞 = 𝑞𝐻(𝑞)− 𝐾(𝑞) 41 where H(1)=ζ(1) is the Hurst exponent for the monofractal case. For monofractals, that is when C1=0 or α=0, ζ(q) is a linear function of q. In the present study, H, C1 and α were estimated in Step 3 by the nonlinear regression of equation (2.7) (e.g Seuront et al. 2002, Tennekoon et al. 2005). 2.2.4 Estimation of confidence intervals We used a nonparametric bootstrap procedure in which the residuals obtained from the fit of the nonlinear regression to the Universal Multifractal Model (equation 2.8) were randomized (e.g. Crawley 2007). The 95% confidence intervals (Efron and Tibshirani 1993) were estimated using 1000 resamples. Calculations were performed using the R software (R Development Core Team 2005, version 2.7.2). 2.2.5 Detection of scaling ranges Scaling ranges were detected on Step 3 using three different criteria: local slopes, zero-slope and “R2-SSR” (Takalo et al. 1995, Seuront et al. 2004). First, we used an automatic adaptive scaling range selection algorithm (Xia et al. 2005) which finds the linear region of a curve by comparing how closely a linear regression model characterizes the observed data points. It uses an independent measure of goodness of fit known as Theil’s inequality coefficient U (Theil 1972), which is calculated as: 2.9 𝑈 = 1𝑁∗ (𝑌! − 𝑌!)!!∗!!!1𝑁∗ (𝑌!!∗!!! )! + 1𝑁∗ (𝑌!)!!∗!!! 42 where Ŷr are the values fitted by the linear model, 𝑌r correspond to the structure function values and N* is the sample size for U equivalent to N/4. The observed time scale or time lag r was measured with increments of one fishing day. The denominator of equation (2.9) is used to normalize U between 0 and 1. As U approaches zero, the fitted data points get closer to the observed data. To determine scaling breaks, we used a threshold value of U equal to 0.015 as estimated between two time scales r1 and r2 (Xia et al. 2005). The initial starting point for the algorithm was set as the middle region of the curve; i.e., r1 = (N*/2)-1 and r2 = (N*/2)+1. Second, to check the consistency of the scaling range selected with the Xia et al. (2005) algorithm, local slopes of Yr were calculated between adjacent data points and plotted as a function of time lag r on a logarithmic scale. If the time series is self-similar, local slopes should be approximately constant for some time period (Takalo et al. 1995) which is equivalent to the so called “zero-slope” criterion used by Seuront et al. (2004). This criterion was designed to detect a difference between the slope of a regression line and a theoretically expected slope of zero using a t-test (Zar 1999). Finally we applied the “R2-SSR” criterion (Seuront et al. 2004), which determines the scaling range maximizing the coefficient of determination (R2) and minimizing the total sum of squared residuals (SSR) of the regression line. 2.2.6 Selection between monofractal and multifractal models To determine whether structure function observations calculated using equation (2.4) are better fitted by a monofractal than a multifractal model, or vice versa, visual detection is generally used. This is done by visually comparing the fits of the theoretically fitted 43 monofractal curve (straight line) and the fitted multifractal curve (dome shaped curve) to the observations. However, we used the likelihood ratio test for nested models (Burnham and Anderson 2002) to select the model which best describes the structure function observations. By setting C1=0 or α=0 in equation (2.7) for the multifractal model, we obtained the monofractal model ζ(q)=qH, denoted as M1, which is a nested model of the multifractal model M2. Assuming normal errors with constant variance the statistic R’ (Hilborn and Mangel 1997), which represents twice the difference in negative log-likelihoods given the data, also called the deviance, was calculated as: 2.10 𝑅! = 2 log𝑳(M! data − log𝑳(M! data)) where L(Mk | data) is the likelihood of model k. R’ has a χ2 distribution with p2-p1 degrees of freedom, where p1 and p2 are the number of parameters of M1 and M2, respectively; here, p2-p1=2. If R’ is greater than χ2 at the α’=0.05 level, then the multifractal model is significantly better than the monofractal model at this level. 2.3 Results 2.3.1 Scaling of TDC and stationarity of TDCr time series A preliminary indication of scaling over multiple ranges was obtained by fitting a straight line through the power spectrum plot of the linearly detrended TDC time series (Fig. 2.4). However, there was significant variability in this plot, especially at both ends of the spectrum, which means that it was not possible to detect clear scaling breaks that can be used to estimate monofractal parameters. This power spectrum analysis is only used here as a preliminary assessment of scaling in the TDC time series (Step 1, Fig. 2.3) and to classify the TDCr time series as stationary or non-stationary (Step 2, Fig. 2.3). 44 Figure 2.4. Power spectrum of the pink shrimp (Pandalus jordani) total daily catch time series (TDC) off the west coast of Vancouver Island for 1994-1996. The fitted straight line is a preliminary indication of scaling. The slope (-β) obtained from the power spectrum analysis (equation 2.3) is equal to -0.37, which means that the TDCr time series corresponds to a stationary process (Davis et al. 1994). A more detailed analysis using structure functions revealed the existence of two increasing trends in TDCr time series; the first between 1 and 3 fishing days and the second between 4 and 15 fishing days (Fig. 2.5). These scaling breaks are in agreement with those observed in the power spectrum plot of the TDC time series at the same time lags (Fig. 2.4). Therefore, the nonstationary segment over the range from 1 to 15 fishing days was excluded before the estimation of multifractal parameters (step 2 of figure 2.3). This latter procedure avoids a biased estimation of the Hurst coefficient due the spurious effect of trends on the self-similar properties of total daily catches. This means that 16 fishing days was selected as the lower limit before the detection of scaling ranges. For the TDCr time 45 series we obtained a value of H=ζ(1)=0.05, after excluding the initial non-stationary segment, which confirms: (i) that TDCr is close to a strict stationary process for which H=ζ(1)=0, and (ii) the need to take cumulative sums of (integrate) the TDCr time series before calculating the structure functions. As explained in the Methods section, taking cumulative sums is necessary when a time series is analyzed as a random walk and H=ζ(1) is close to zero (Gao et al. 2006). The integrated TDCr time series forms the total daily catch anomaly time series (TDCa) and the latter is used hereafter in the calculation of scaling ranges and in the estimation of multifractal parameters. Figure 2.5. Structure function Yr(q) versus r in a log-log plot obtained from the pink shrimp (Pandalus jordani) total daily catch residual time series (TDCr) for the west coast of Vancouver Island from 1994-1996. Two increasing trends, the first between 1 and 3 fishing days and the second between 4 and fishing 15 days are visible. 46 2.3.2 Selection of scaling ranges The Xia et al. (2005) adaptive search algorithm applied to detect the linear region of the curve using structure function values Yr(q) for q=1 selected a scaling range of 16 to 120 fishing days (Fig. 2.6). The value of Theil’s inequality coefficient U obtained from equation (2.8) for the selected scaling range was 0.0001, which is far below the threshold value of U=0.015. A value of R2=0.99 (p<0.01) confirmed the good fit of the linear regression model for this scaling range. Figure 2.6. Scaling range (between 16 and fishing 120 days) selected by the adaptive search algorithm (Xia et al. 2003) applied to detect the linear region of the structure function curve of the total daily catch anomaly time series (TDCa) for q=1. Upper scaling limit is marked by an arrow. Inset shows a close-up view of the departure of the structure function from the linear trend at 120 fishing days. According to the second criterion for detecting scaling ranges, local slopes were plotted against the time lag and a zero slope was obtained for several regression lines. These lines have a fixed lower scaling limit of 16 fishing days and an upper scaling limit that varied 47 between 116 and 129 fishing days. After approximately 120 fishing days, local slopes generally decreased showing the absence of a scaling region after that time lag (Fig. 2.7). The response of structure functions to the presence of periodic oscillations is well known. In these cases, structure functions reach their maximum values at half of their period and their minimum values at the time scale of their period (Takalo et al. 1995, Yu et al. 2003). Figure 2.7. Temporal variability of local slopes calculated for TDC anomaly time series (TDCa) for time lags between 16 and 321 fishing days. Local slopes fluctuate around a constant value up to approximately 120 fishing days and notably decrease after that. After stitching together the data points of the total daily catch time series from the main fishing seasons from 1994 to 1996, a period equal to 214 fishing days was artificially created. The structure function for the TDCr time series reaches a value near its maximum at 214 fishing days which demonstrates that the main seasonal period of 214 fishing days 48 was efficiently removed and has no effect on the structure function plot. Using the “R2-SSR” criterion and a lower scaling limit of 16 fishing days we obtained a maximum coefficient of determination (R2=0.9978) and a minimum sum of squared residuals (SSR=0.0475) for a regression line with zero slope (t=0.621, p<0.01) at a time lag equal to 120 fishing days. In conclusion, lower and upper scaling limits of 16 and 120 fishing days were used to estimate monofractal and multifractal parameters. 2.3.3 Estimation of monofractal and multifractal parameters and selection of models Structure functions for the TDCa time series exhibit linear trends over the selected scaling range (16 to 120 fishing days according to Xia’s et al. (2003) algorithm) for all statistical moments. For clarity, only four moments (q=1, 2, 3, 4) from a total of 40 are shown in Fig. 2.8. The value estimated for H was 0.8779 (95% CI: 0.8777, 0.8783). This value can also be obtained using equation (2.6) for which H=ζ(1)=0.8778 which is very close to the value from the nonlinear regression model. The empirical structure function calculated using the exponents ζ(q) for q between 0.1 and 4.0 is non-linear which means that the TDCa time series is multifractal over the selected scaling ranges (Fig. 2.9). The empirical and theoretical multifractal curves also showed good agreement (Fig. 2.9) when the theoretical curve was calculated using the estimated universal multifractal parameters C1=0.0429 (95% CI: 0.0425, 0.0435) and α=1.5073 (95% CI: 1.493, 1.521). Since the TDCa time series is multifractal, H=ζ(1) cannot be interpreted as a Hurst coefficient because multiple scaling patterns are required to characterize the dynamics of this time series. This means that a single coefficient (H) cannot determine the complete statistical scaling features of the TDCa time series, and additional parameters (C1 and α) are necessary to explain the 49 variability. For this case, the Hurst coefficient is estimated using the Lévy index as H=1/α (Samorodnitsky and Taqqu 1994) which is equal to 0.66. Figure 2.8. Structure function Yr(q) versus r in a log-log plot obtained from the pink shrimp (Pandalus jordani) total daily catch anomaly time series (TDCa) for the west coast of Vancouver Island from 1994-1996. A linear trend selected by the adaptive search algorithm (Xia et al. 2003) for q=1 in Figure 2.6 is visible also for statistical moments q=1, 2, 3, 4 for a range of scales from 16 to 120 fishing days marked between vertical dotted lines. The estimated R’ statistic from the likelihood ratio test obtained from equation (2.10) was 362; this is much greater than 5.99, the value calculated at the α’=0.05 significance level. This test clearly selected the multifractal over the monofractal model, which confirms the good fit of the theoretical multifractal model to the empirical multifractal observations shown in Fig. 2.9. However, more observations will be necessary to make multifractal models a useful tool in the forecasting of total daily catches in fisheries research. This is 50 because the error introduced in using finite length time series for the estimation of scaling exponents using structure functions is less than 5% if the minimum number of data points varies between 1,000 and 3,000 according to the dynamics of the time series (Kiyani et al. 2009). This does not mean that the TDCa time series is not multifractal. This only means that our estimation of multifractal parameters could have been affected by the limited number of observations available, and in consequence when using structure functions and a multifractal model to forecast TDC values no less than 3,000 observations should be used. Figure 2.9 Empirical curve for scaling exponents, ζ(q), (filled circles) compared to the monofractal curve ζ(q)=qH (continuous line) and to the theoretical universal multifractal function (crosses) obtained from the pink shrimp (Pandalus jordani) total daily catch anomaly (TDCa) time series with H = 0.8779, C1 = 0.0429 and α =1.5073. 2.4 Discussion 2.4.1 Identification of multifractal patterns To our knowledge, this is the first time that multiscale temporal patterns of variability, i.e. multifractality, have been demonstrated for an invertebrate fishery. Previous studies 51 conducted on invertebrate fisheries found monofractal spatial patterns, e.g. fractality in the spatial variability of catch rates and abundance estimates from the Northern Prawn Fishery in Australia (Wang et al. 1999) and from the krill fishery in the Southern Ocean (McClatchie et al. 1994). Monofractal patterns have also been found in interannual growth rates and recruitment time series of several fish and invertebrate stocks (Niwa 2006, 2007). However, these studies focused only on the estimation of the Hurst parameter to quantify the long-range correlation structure in these time series, and did not consider the possible existence of multiple scaling patterns that could produce more complex dynamics. Our findings support those reported by Montes (2004) who analyzed the temporal multifractal dynamics of catches and catch per unit effort from the southern hake (Merluccius australis) fishery in southern Chile. This pattern of self-similar temporal dynamics may be a widespread feature of fisheries, at least in the Eastern Pacific Ocean, and not restricted to specific ecosystems or particular fisheries. Our results demonstrate the existence of multifractal dynamics in smooth pink shrimp daily catch time series off Vancouver Island, Canada, with a scaling range of 16-120 fishing days. This scaling range suggests that different processes affect and structure the availability/catchability of P. jordani at time scales less than 120 fishing days, compared with longer time scales. As different segments of the TDCa time series have different scaling properties, the simplest explanation is that low and high catches are regulated by different mechanisms, rather than a single mechanism. Several mechanisms have been proposed to explain the emergence of fractal-like spectra in the abundance variability of terrestrial and aquatic populations, including fish stocks, at an inter-annual scale. For example, age structure, density-dependence, measurement errors and the interaction 52 between age structure and stochastic recruitment are factors used to explain the predominance of low-frequency over high-frequency variability (e.g. Akçakaya et al. 2003, Bjørnstadt et al. 2004, Halley and Inchausti 2004) that lead to the appearance of long-range correlations and the characteristic linear trend given by equation (2.3) of the power spectrum plot. However, according to our results, intra-seasonal time scales of variability play major roles in structuring the scaling dynamics in this pink shrimp fishery. The fractal approach recognizes intermediate scales of variability in this fishery (i.e. from weeks to a season) as fundamental in the emergence of multiple scaling patterns in pink shrimp catches. Note that only fishing days were analyzed over the main fishing seasons from 1994 to 1996. Since TDCr was demonstrated to be a stationary time series across these three fishing seasons, the patterns of variability in any one season are consistent with those in the other seasons, despite the intervening November to March calendar (and non-fishing) days. As a consequence, the interpretation of the processes underlying the 16 to 120 fishing days scaling range apply to the spring to fall fishing seasons. 2.4.2 Potential applications of fractal models in the pink shrimp fishery For a multifractal time series, the scaling symmetry can be used to predict the probabilities of extreme events at time scales beyond the longest available or at magnitudes in excess of those in the current observations (e.g. Finn et al. 2001). In this latter case, the scaling of the time series involves the preliminary estimation of H, α and C1. For example, using the multiple scaling patterns observed for the pink shrimp total daily catch time series it is possible to predict the probability of achieving values two, three or other multiples of the 53 observed mean for this time series (e.g. Schertzer et al. 1997, Finn et al. 2001). In terms of its application to the management of this fishery, such predictions of extreme variability can be used to assist fisheries managers with explaining and interpreting occasional very large catches. More detailed information about the multifractal theory used for the prediction of extreme events can be found for example in Schertzer et al. (1997) and Finn et al. (2001). Up until 1997, the pink shrimp fishery along the west coast of Vancouver Island was managed primarily by seasonal closures. Since 2000, management approaches have included preliminary catch ceilings based on pre-season biomass forecasts and the application of pre-determined harvest rates, followed by re-assessment and revision of Total Allowable Catches after completion of fishery-independent biomass surveys (Rutherford et al. 2004; DFO 2011). When the Total Allowable Catch for a particular area is reached, that area is closed. It is not known at the start of a season whether, and when, a particular area may reach its quota and be closed. The universal multifractal parameters H, α, and C1 can be used to construct fractal algorithms (Schertzer et al. 1997) to forecast total daily catch values at intra-annual time scales. These algorithms take into account the extreme fluctuations and the scaling symmetry of multifractal time series and can significantly reduce the forecast errors of individual data points (Richards 2004). In addition, the multiple scaling patterns for this time series can potentially provide an in-season estimate of the future variability of catches using historical data and observations from the early part of the season. Taken together, these properties could provide an early season estimate of when the catch ceiling may be reached in an area, and the variability (i.e. uncertainty) around this date (assuming that the fishery and environmental conditions 54 proceed as they have in the past). Traditionally, short-range memory models (ARIMA) have been used in fisheries research to obtain forecast values of catches (Ryall et al. 1999, Hanson et al. 2006, Tsitsika et al. 2007) or a long-memory parameter is used to analyze the dynamics of catches (Garcia et al. 2010), but these fail to represent the extreme variability of catches on a daily scale. In addition, if the observed patterns of variability differ from the predicted patterns, it may indicate that something has changed in the underlying processes governing daily catches, thereby serving as an early-warning indicator of these changing conditions. By using fractal and multifractal models it is also possible to detect the degree of association between fisheries and environmental/oceanographic time series in a more formal way than using classical correlation analyses. This new correlation function based on multifractal theory has desirable properties: (i) no assumptions about the probability distribution of each time series are made, (ii) the range of scales for which a significant correlation is found is detected, taking into account the scaling symmetry of each time series, and not imposed (as usually done) prior to the analysis, (iii) the degree of association is calculated considering different intensity levels (by using different statistical moments q) and not only the mean of each time series (Seuront and Schmitt 2005a). This latter property can be used to track the evolution of the correlation between two time series. For example, in a hypothetical fish population high catches can be positively correlated with high and intermediate levels of sea surface temperature (SST), but not associated at all or negatively correlated with low values of SST. Detailed examples of this application in oceanographic research can be found in Seuront and Schmitt (2005a,b). Preliminary analyses conducted using much longer time series of catches (N~2500 observations) confirm the existence of 55 long-range cross-correlations between pink shrimp TDC and sea surface temperature (SST) (not shown here). 2.4.3 Hypotheses for the causes of these fractal patterns Seasonal scaling patterns may arise from biological processes regulating the abundance and availability of shrimp to these fishing grounds, from behavioural processes within the shrimp fishery itself, from management actions, or a combination of all. Having demonstrated the existence of multifractal temporal patterns in this fishery, ongoing work is focused on quantifying the effects of physical factors and management actions on the dynamics of catches from a fractal perspective. Various hypotheses are being tested and the results will be presented in forthcoming publications. In the following paragraphs we identify some of these hypotheses which may account for the upper scaling limit and multifractality in this fishery. For biological processes, the upper scaling limit (i.e. 120 fishing days) of the TDCa time series may be associated with the characteristic time scale of seasonal migrations of important predators of these shrimp. Pacific hake (Merluccius productus) coastal populations, one of the predators of P. jordani (Buckley and Livingston 1997), undertake a south-north seasonal migration through the coastal upwelling domain of the North-East Pacific for feeding purposes (Beamish et al. 2005). It is not clear how long hake remain in Canadian waters off the west coast of Vancouver Island but they are thought to be present between 80 and 180 days during warm years (Martell 2002). Food web processes relating to the influence of meanders and eddies on phytoplankton production with seasonal time scales of 100-200 days (Henson and Thomas 2007) could also influence the abundance and 56 availability of pink shrimp to the fishery. These values encompass the upper scaling limit (120 fishing days) found for this pink shrimp catch time series, although more evidence is necessary to support these hypotheses. In regards to processes within the fishery itself the combined effect of multiple factors like an increase in the abundance of shrimp (Rutherford et al. 2004), and an increase of fishing effort (DFO 2001) may explain the emergence of multifractal scaling patterns. Significant changes in the level of multifractality in time periods with different fleet-resource dynamics were detected in the southern hake fishery of Chile (Montes 2004). To explore this hypothesis, analyses should be conducted which incorporate, for example, time periods for which the shrimp population and fishing effort were more homogeneously distributed compared with those periods when significant declines in abundance occurred. Socio-economic factors can have an effect on the selection of fishing grounds and on the behaviour of a trawl fleet (Walters and Martell 2004), which can cause short-range correlations. For example, in the case of rose shrimp (Aristeus antennatus) in the NW Mediterranean the better knowledge of shrimp distributions obtained by the fleet in the course of the week, and higher prices obtained by the fishermen at the end of each week, were proposed to explain significant short-range autocorrelations in catch rates at a time lag of 5 days (Sardà and Maynou 1998). In our study, scaling breaks were detected for pink shrimp total daily catch time series at 3 fishing days which may reflect the effect of fishing location and the duration of the fishing trips. For the P. jordani fishery off the coast of California, the locations of fishing on any given day were determined in part by high catches the previous day (Eales and Wilen 1986). In this latter study, fishers used this information to choose the same fishing location until the high density of shrimp was 57 reduced. In the case of the Pacific hake fishery off the west coast of North America, fishing persisted at specific location for 3-4 days (Dorn 2001). These examples suggest that socio-economic factors may have an effect on the variability of catches but restricted to time scales no greater than a week. For the P. jordani fishery off the west coast of Vancouver Island, wind stress and tidal speed, among other environmental factors, can affect the availability of this species on a time scale of 8-10 days (Perry et al. 2000). This scale is somewhat shorter than the range obtained in this study (16 to 120 fishing days), however, this does not mean that environmental variability has no effect on the availability on P. jordani on longer time scales. In the study conducted by Perry et al. (2000), using data from 1996, the variability in catch rates of P. jordani was analyzed up to a maximum of 90 days, and indeed, low frequency oscillations in catch rates were found at 30-40 days. Our ultimate goal is to identify underlying processes that may be sources of the observed fractal properties, as a guide to understanding catch fluctuations and to improve the management of this fishery. These scaling patterns could be used to construct early warning indicators of change in this, and possibly other, fisheries. However, fractality needs to be demonstrated for other fisheries operating in the same and other ecosystems in order to determine its generality and to be accepted as an early warning fisheries indicator. Considering the increasing monitoring activities of fisheries operations, including on-board observers, the fine resolution data necessary to detect and quantify scaling patterns should become more readily available. Here, we have demonstrated the existence of multiple scaling patterns in pink shrimp catch time series off the west coast of Vancouver Island, which is the necessary first step to demonstrate the ubiquity of fractal patterns in fisheries 58 time series. Having identified such properties in this study, forthcoming papers will analyze longer time series of smooth pink shrimp catches including periods with different management actions and abundance levels, to examine the use of these properties as early indicators of changes in fisheries. There is an urgent need for indicators that anticipate changes in populations of exploited marine species. Complex patterns of variability found in time series that are not detected using classical models in fisheries were modeled. Fractal analyses were used to detect time-invariant scaling symmetries in daily catch time series from the smooth pink shrimp (Pandalus jordani) fishery from the west coast of Vancouver Island, Canada. A universal multifractal model, which accounts for intermittent fluctuations and extreme values in a time series, provided a better fit to daily catches than a monofractal model. Multifractality is an indication of multiple scaling patterns and suggests that more than one process is affecting the variability of catches. Fractal dynamics in catch time series were found for a range of scales between 16 and 120 fishing days. To our knowledge, this is the first time that multifractality has been demonstrated for an invertebrate fishery. The multifractal model has the potential to provide an in-season estimate (up to 120 fishing days) of the variability of shrimp catches based on the variability at time scales less than one month. Changes in these fractal patterns may provide an early warning that conditions underlying this fishery are changing. In conclusion, smooth pink shrimp total daily catches (TDC) follow a multifractal temporal dynamics during 1994-1996 fishing seasons (N=642 fishing days, 3 stitched fishing seasons). In order to see if this holds true for the longer time series, in Chapter 3 we analysed catch dynamics between 1987 and 1998 (N=2568 fishing days, 12 stitched fishing 59 seasons). Unlike this Chapter, in Chapter 3 the Detrended fluctuation analysis (DFA) methodology was used because it allowed the estimation of fractal parameters in a time series with missing observations. Overall, 16.5% of daily data were missing observations. Independently, fractality and multifractality were analyzed in oceanographic time series, including sea surface temperature and wind stress. At the end of Chapter 3, bivariate fractal analyses were performed to quantify how the fractal dynamics of sea surface temperature (and wind stress) were affecting the temporal fractal pattern of variability in smooth pink shrimp catches. 60 Chapter 3: Estimation of long-range temporal dependence of sea surface temperature and wind stress and their effects on smooth pink shrimp (Pandalus jordani) fisheries catches off the west coast of Vancouver Island 3.1 Introduction Long-range dependence (LRD), also called long-range persistence (LRP) or long-range memory, is present in numerous hydro-climatic and geophysical time series (Pelletier and Turcotte 1999, Kantelhardt et al. 2006, Koscielny-Bunde et al. 2006, Ghil et al. 2011, Rybski et al. 2011, Lovejoy and Schertzer 2013, Witt and Malamud 2013). Furthermore, LRD has been incorporated into climate variability models (Fraedrich and Blender 2003, Vyushin and Kushner 2009) and also used as a statistical tool to monitor the climate (Kurnaz 2004, Mudelsee 2010). In oceanographic research, LRD has not received much attention until the last decade when long-range memory (fractionally differenced models) and short-range memory (autoregressive integrated moving average, ARIMA models) were compared quantitatively as to their ability to explain oceanographic regime shifts and inter-decadal variability in the North Pacific Ocean (Percival et al. 2001, Rudnick and Davies 2003, Overland et al. 2006, Khaliq and Gachon 2010). During the last decade, several oceanographic studies have detected and quantified LRD in time series of upwelling (Denny et al. 2004), ocean surface velocities (Isern-Fontanet et al. 2007), chlorophyll (Zhan 2008) and sea level (Barbosa et al. 2008, Dagendorf et al. 2014). Long-range dependence (LRD) for sea surface temperature time series was detected in the North Pacific Ocean (Zhu et al. 2010), off the California coast (Lewis and Ray 1997), on east and west sides of the tropical Atlantic ocean (Royer 61 and Fromentin 2007), in the North Atlantic Ocean (Monetti et al. 2003, Rouyer et al. 2010, Seuront 2010) and off the west coast of Africa (Mendelssohn 1990, Nieves et al. 2007). The existence of LRD in time series of sea surface wind stress and/or ocean transport is less documented but some evidence is presented by Mendelssohn (1989), Yano et al. (2004) and Ashkenazy and Gildor (2009). Usually, oceanographic studies focused on the estimation of long-range dependence parameters, but recently, LRD has also been used for the simulation of oceanographic responses to atmospheric forcing (Legatt et al. 2012). Long-range dependence in the dynamics of fish populations is frequently associated with changes in low-frequency oscillations of oceanographic variables. This means that changes in the temporal dynamics of fish catches, catch per-unit effort or spawning stock biomass can arise due to the emergence of oceanographic low-frequency components (Rouyer et al. 2010, Lindegreen et al. 2013). In these latter studies, oceanographic low-frequency oscillations were considered as the main source of long-range dependence (LRD), but it is important to remark that from a statistical standpoint LRD can emerge without forcing by oceanographic conditions. Smooth pink shrimp fishery management areas off the west coast off Vancouver Island, Canada, are imposed on a complex and dynamic oceanographic transition zone where the North Pacific Current (NPC) diverges into the northward flowing Alaska Coastal Current (ACC) and the southward flowing California Current (CC) (Thomson 1981, Perry et al. 2007). This transitional zone is strongly influenced by climate cycles superimposed one on another, whose interaction is not completely understood (Patterson et al. 2013). Among these cycles, one of the most recognized is the El Niño phase of the El Niño Southern 62 Oscillation (ENSO), which directly affects water properties altering the distribution of species in this ecosystem (Okey et al. 2014). The spring to fall period each year, which also corresponds to the smooth pink shrimp (Pandalus jordani) fishing season off the west coast of Vancouver Island (Perry et al. 2000), is characterized by predominately seaward wind-forced Ekman transport of the surface layers and upwelling conditions due the southward displacement of the North Pacific high pressure system (Hickey 1998). Specific features of upwelling off the west coast of Vancouver Island are mainly driven by winds and the complex local topography (Fang and Hsieh 1993). Long-range dependence in oceanographic time series can be caused by non-stationary low-frequency scales of variability, such as the El Niño Southern Oscillation (ENSO, Boucharel et al. 2009) or the Pacific Decadal Oscillations (PDO, Khaliq and Gachon 2010). For the North-East Pacific ecosystem, significant shifts in climate indices representing sea surface temperature variability are associated with long-range memory process for which autocorrelation extends over multiple years (Overland et al. 2008, Yuan et al. 2014) causing changes in the temporal dynamics of fish populations (Litzow and Mueter 2014). For this ecosystem, three well resolved dominant low frequency scales of variability have been detected: quasi-biennial oscillations (QBO, 2-3 years), El Niño Southern Oscillation (ENSO, 5-7 years) and Pacific Decadal Oscillations (PDO) with phases persisting for 20-30 years (Ware 1995, Hollowed et al. 2001, Bailey et al. 2004, Möllmann and Diekmann 2012). All of these can contribute to the emergence of long-range dependence in time series of fish populations. Under this framework, the first goal of this Chapter was to detect statistical long-range dependence (LRD) for sea surface temperature and wind stress from the west coast off 63 Vancouver Island. The second goal was the detection of long-range dependence in smooth pink shrimp total daily catches, and to verify if its scaling range -previously detected in Chapter 2- varies using a longer time series of catches. To accomplish this goal, a new method known as Detrended Fluctuation Analysis (DFA) was applied. The third objective was to detect if the fractal dynamics of Pandalus jordani catches is affected by the fractal dynamics of oceanographic time series. Fisheries and oceanographic time series are characterized by several processes acting on different time scales (Rouyer et al. 2008), which make the detection and quantification of long-range dependence a difficult task. For that reason, a well-known method for the estimation of long-range dependence was selected: detrended fluctuation analysis (DFA; Peng et al. 1992, Kantelhardt 2011, Witt and Malamud 2013), due to its ability to estimate LRD from complex non-stationary time series (Chen et al. 2002, Kantelhardt 2009). In addition, DFA allows the estimation of fractal parameters in time series contaminated with data gaps. 3.2 Methods 3.2.1. Oceanographic databases Daily values of sea surface temperature (SST), measured at Amphitrite Point Lighthouse off the west coast of Vancouver Island (48º55’N; 125º32’W; Fig. 3.1) as part of an observational program covering shore stations along the British Columbia coast (Hollister and Sandnes 1972), were obtained from Fisheries and Oceans Canada (http://www.pac.dfo-mpo.gc.ca/science/oceans/data-donnees/lighthouses-phares/index-eng.html). Sea surface temperature records were obtained within an hour of local daytime high tide from the top 64 meter of the water column (Cummins and Masson 2014). Figure 3.1. Amphitrite Point Lighthouse (48º55’N; 125º32’W, circle) and buoy number 206 (48.84 ºN; 126º00’W, triangle) located off the west coast of Vancouver Island, Canada. This database covers the period 1935-2013 (79 years) at a daily scale, but only observations from 1940 to 2013 (N=27,010 daily observations; 74 years) were analyzed. Leap days were omitted, and a total of 969 missing observations (approximately 3.6%) were detected (Fig. 3.2), which were replaced by a long-term mean calculated for that specific day (Talkner and Weber 2000). For example, if no observation was registered for April 1, 1987, the average value calculated considering every 1st of April from 1940 to 2013 was used to fill this missing observation. We are aware that it is also possible to estimate long-range dependence parameters for time series with missing values without the need to fill missing 65 observations (Chamoli et al. 2007). However, this implies the loss of temporal sequence between observations and for that reason, missing values were estimated to avoid mis-interpretations of results from bivariate fractal analyses. Figure 3.2. Sea surface temperature time series (SST in °C) at Amphitrite Point Lighthouse for the 1940-2013 period. Missing values (3.6%) were replaced by a long-term mean (see text) and their positions marked by red crosses centered at 5ºC. A horizontal line is centered at time series mean (10.10ºC). Vertical lines separate periods of ten years. The daily marine surface wind database, calculated from meteorological buoy number 206 located 28.4 km off the WCVI (48.84 ºN; 126.00ºW, Fig. 3.1) by Faucher et al. (1999), was obtained from R.I. Perry (personal communication). This database corresponds to the wind stress (i.e., square of the wind speed), which extends from 1958 to 1998 (41 years) at 6 hourly intervals with no missing observations. Faucher’s wind reconstructions can be considered a reliable wind database, and for that reason, it can be used to calculate an 66 accurate index of coastal upwelling (Foreman et al. 2011). The wind database was aggregated on a daily scale, which contains N=14965 daily mean values of wind stress (WS) (Fig. 3.3a). Figure 3.3. a) Time series of daily mean wind stress (WS, m/s2) for the 1958-1998 period and its mean annual cycle (green line). b) Residual wind stress time series. Vertical lines separate periods of ten years. 3.2.2 Fishery database Total daily smooth pink shrimp catch time series (TDC) for 1987-1998 stitched fishing seasons (see Methods, Chapter 2) (N=2568, 12 years, 16.5% of missing observations) were calculated (Fig. 3.4a). 67 Figure 3.4. a) Time series of smooth pink shrimp total daily catches (TDC, 1000 Kg) for the 1987-1998 period, and b) TDC anomaly time series. Missing values (16.5%) were replaced by a long-term mean (see text) and their positions marked by red crosses centered at 1000 Kg. Vertical grey lines separate periods of 428 days, which correspond to two continuous fishing seasons. No fishing activity for a specific day of the fishing season means that a missing value for that day was obtained. In order to be consistent with the filling of missing observations applied to oceanographic time series, the same procedure was applied to the TDC time series; i.e., missing TDC observations were replaced with its long-term mean calculated for each day of the 12 fishing seasons. Definitions of variables used in this chapter are summarized in Table 3.1. 68 Table 3.1. Definitions of variables and parameters used in Chapter 3. Symbol Definitions LRD Long-range dependence LRP Long-range persistence ARIMA Autoregressive integrated moving-average model AR-1 Autoregressive process of order 1 NPC North Pacific Current ACC Alaska Coastal Current CC California Current ENSO El Niño Southern Oscillation PDO Pacific Decadal Oscillation PNA Pacific North American pattern SST Sea surface temperature QBO Quasi-Biennial Oscillation DFA Detrended Fluctuation Analysis SR Scaling range 69 Symbol Definitions N,n Number of observations WCVI West Coast off Vancouver Island TDC Total daily catches TDCa Total daily catch anomalies WS Wind stress CWT Continuous wavelet transform x’n Reconstructed time series δj, δt Scale and time intervals Cδ Reconstruction factor ℛ{Wn(sj)} Real part of continuous wavelet transform sj Wavelet scales j1, j2 Lower and upper wavelet scales SSTF Filtered sea surface temperature WSF Filtered wind stress Y(i) Profile function (cumulative sum) 70 Symbol Definitions N/s Number of segments s Length of segment v Defined segment pv Polynomial function fitted to segment v m Order of polynomial function σ Variance F(s) Fluctuation function for segment of length s DFA-1 Detrended fluctuation analysis of order 1 α Slope estimated from fluctuation function plot 1/f Flicker noise (α=1) MF-DFA Multifractal Detrended Fluctuation Analysis q Statistical moment H Hurst exponent h(q) Generalized Hurst exponent o,f Oceanographic and fishery time series 71 Symbol Definitions Fo,f Bivariate fluctuation function Fo,f(q,s) q-order detrended covariance ho,f Generalized Hurst exponent for bivariate multifractals β Power spectrum slope fGn Fractional Gaussian noise fBm Fractional Brownian motion (summed fGn) HML Maximum likelihood estimates for the Hurst coefficient smax Upper scaling range HDFA Detrended fluctuation analysis estimates of Hurst coefficient 𝐻!" Average value of Hurst coefficients L(x,β) Maximum likelihood function maximized for time series x |R(β)| Determinant of the covariance matrix defined for exponent β xT Transposed time series x 72 3.2.3 Removal of unwanted oscillations Periodic oscillations, considered as deterministic trends which distort the estimation of long-range dependence (LRD) parameters, were removed (Montanari et al. 1999, Gencay et al. 2002). For sea surface temperature (SST) and wind stress (WS) time series the mean annual cycle (365 days period) was subtracted from original observations forming SST and WS anomalies, respectively (Kiss et al. 2007, Franzke 2012). The mean annual cycle for the smooth pink shrimp time series, which has a length of 214 days formed by the stitching of continuous fishing seasons (see Chapter 2) was removed, and the new time series defined as TDC anomaly (Fig. 3.4b). The filtering of trends of sea surface temperature and wind stress anomalies to isolate their stochastic variability and long-range dependence properties was also necessary (Yuan et al. 2014). Daily sea surface temperature anomalies, calculated using N=74 years (1940-2013), were first normalized (mean removed and divided by standard deviation, Dakos et al. 2012), and then periodicities between 128 and 4096 days (0.35 and 11.22 years, respectively) were removed. These latter lower and upper limits were selected because this range includes the sea surface temperature annual cycle or first harmonic (365 days) and its second harmonic of 182 days (Gower 2002). The removed frequency band also contains well-known periodic or pseudo-periodic SST oscillations in the North East Pacific, such as the quasi-biennial oscillations (2-3 years) and El Niño Southern Oscillation (5-7 years) (Ware 1995). These periodicities were reconstructed and filtered out using a continuous wavelet transform (CWT, Torrence and Compo 1998) according to: 73 3.1 𝑥!! = 𝛿!𝛿!!/!𝐶!𝜓!(0) ℛ 𝑊! 𝑠!𝑠!!/!!!!! where 𝑥!! correspond to the reconstructed time series, 𝛿! and 𝛿! are the scale and time intervals, 𝐶! is a reconstruction factor equal to 0.776 when the Morlet mother wavelet (𝜓!) is used. The term ℛ{Wn(sj)} is the real part of the wavelet transform defined over scales sj, and j1 and j2 correspond to selected lower (128 days) and upper (4096 days) scales, respectively. The same procedure was applied to daily mean wind stress time series (WS), which means that individually reconstructed time series 𝑥!! were subtracted from SST and WS anomaly time series, respectively. After filtering out periodic oscillations, they are called from here on as filtered sea surface temperature (SSTF) and filtered wind stress (WSF) time series. It is necessary to remark that when oceanographic time series were filtered with the objective to test their bivariate fractal relationship with the smooth pink shrimp time series, only daily observations matching fishing seasons (April-October) from 1987 to 1998 were used. This means that as N=2568 daily values were used to filtered out periodic components, the upper selected wavelet scale j2 defined in the above paragraph decreases from 4096 days (11.22 years) to 2048 days (7.03 years). 3.2.4. Fractal analysis Long-range dependence (LRD) was detected and quantified using Detrended Fluctuation Analysis (DFA) methodology, used and recommended to estimate long-range dependence in complex geophysical time series (Kantelhardt et al. 2001, Oswiecimka et al. 2006, Shao et al. 2012). DFA has a negligible bias when N>2,000 data points (Delignieres et al. 2006, 74 Fattahi et al. 2011, Witt and Malamud 2013) and also is robust to severe data loss for this latter time series length (Ma et al. 2010). Our shortest time series (fisheries catches, N=2,568) exceeds this lower threshold, and consequently, DFA was considered reliable for the analysis of fisheries catches and oceanographic time series. In the first stage of the DFA methodology, the cumulative sum, also called the nonstationary profile function Y(i), was constructed according to: 3.2 𝑌 𝑖 = 𝑋! 𝑖 = 1,… ,𝑁!!!! where X(k) represents sea surface temperature, wind stress or smooth pink shrimp total daily catches of length N. Then, Y(i) obtained using equation (3.2) is divided into N/s non-overlapping segments of equal length s. In the third stage, for each segment v of length s we fitted a polynomial function 𝑝!(𝑖) of order m, which represent the local trend in that segment. Then, for each segment v of size s (which represent a time varying scale), the variance σ was estimated after subtracting the local trend as: 3.3 𝜎 𝑣, 𝑠 = !! 𝑌 𝑣 − 1 𝑠 + 𝑖 − 𝑝!(𝑖) !!!!! In other words, the fluctuation of the time series is determined as the variance upon the estimated local trend. In the final stage, the variances of all segments were averaged and the square root was taken to obtain the mean fluctuation function F(s) as: 3.4 𝐹 𝑠 = !! 𝜎!!!!! (𝑣, 𝑠) !! In this final stage, we average F(s) calculated over several segments of different size s and analyze the relationship between the mean fluctuation function and the segment size. F(s) 75 obtained using equation (3.4) is only defined for s ≥ m+2, where m is the order of the polynomial (Kantelhardt et al. 2009). Local trends can be estimated by least squares fitting of polynomials of different orders. In this research DFA-1 was used, which means that linear polynomials were removed from time series defined in equation (3.2) A power-law relationship between the fluctuation function F(s) and the segment size s is the signature of scaling and long-range dependence, which was estimated according to: 3.5 𝐹 𝑠 ∼ 𝑠! The scaling exponent or LRD parameter α was estimated by plotting the average root-mean-square fluctuation function F(s) calculated using equation (3.4) versus the segment size s on a double logarithmic scale (Kantelhardt et al. 2001). The slope (α) of the line obtained from the least-squares estimation of equation (3.5) represents the strength of long-range dependence. If α=0.5 there is no correlation and the time series is called white noise; i.e., any observation of the time series is completely uncorrelated from any previous value. A α exponent between 0.5 and 1.0 indicates long-range correlations or LRD (Dijkstra 2013). In other words, a large fluctuation observed in the time series (in relation to its mean value) is likely to be followed by a large fluctuation and vice versa. When 0<α<0.5 the time series is classified as anticorrelated; i.e. large and small values are more likely to alternate (Peng et al. 1995). Values of α ≥1.0 are considered a sign of long-range dependence but not of a power-law form. Special values of α are 1.5 and 1.0. The former correspond to Brownian noise which is the integration of white noise and the latter is called 1/f noise and can be interpreted as a compromise between the total unpredictability of white noise (very noisy signal) and the smooth profile of Brownian noise (Peng et al. 1992). 76 Deviations from linearity are frequently found in the double logarithmic plot when applying DFA methodology at small scales s (s ≈ 10; Kantelhardt et al. 2001, 2002). These deviations are intrinsic to the methodology because the scaling behaviour is only approached asymptotically; i.e., when N approaches a large number. Therefore, the use of small time scales in the fitting procedure will lead to an overestimation of α. Taking into account that scaling breaks at small time scales (s=16 days) have been detected for smooth pink shrimp total daily catch time series (Chapter 2), the use of time scales smaller than 16 days was avoided. Scales (or time series segment size s) up to ¼ of the total time series length N were used. When time scales larger than this threshold were used, the number of segments used in the averaging of F(s) obtained with equation (3.4) became very small, and hence, the calculation of F(s) unreliable (Kantelhardt et al. 2002). 3.2.5. Multifractal analysis The DFA methodology allows the estimation of a single scaling exponent to detect fractal (monofractal) dynamics of a time series. It is also possible to generalize the DFA analysis to analyze time series under a multifractal approach. This method is called Multifractal Detrended Fluctuation Analysis (MF-DFA; Kantelhardt et al. 2002, Kantelhardt et al. 2006) and the first stages are the same as those described for the monofractal DFA analysis in previous paragraphs. After the detrending procedure it is possible to obtain the qth order fluctuation function as: 77 3.6 𝐹! 𝑠 = 1𝑁! 𝐹!! !!!!!!! where q (statistical moment) can adopt any real value. Then the relationship between Fq(s) and the time scale s was analyzed. For the special case q=2 this method is equivalent to the standard DFA (Kantelhardt et al. 2006). If the time series under investigation is long-range correlated, Fq(s) will increase for larger values of s in an analogous way to equation (3.5) as: 3.7 𝐹! 𝑠 ∼ 𝑠!(!) where the exponent h(q) will depend on q. For stationary time series, h(2) corresponds to the Hurst exponent and hence h(q) is designated as the generalized Hurst exponent (Kantelhardt et al. 2002). For fractal (monofractal) time series the h(q) exponent is independent of q. For these types of time series the scaling behaviour of the variances is equal for all time scales and q values. For multifractal time series a single scaling exponent like h(2) is not sufficient to completely characterize their dynamics since many subsets of the time series have different scaling behaviours. 3.2.6 Bivariate multifractal analysis In a similar way as Detrended Fluctuation Analysis (DFA, univariate fractal method) was generalized to obtain the Multifractal Detrended Fluctuation Analysis (MF-DFA, univariate multifractal method), it is possible to generalize DFA to analyze associations between two fractal time series (bivariate fractal method, Podobnik et al. 2007, Podobnik and Stanley 2008) and MF-DFA to analyze associations between two multifractal time series (bivariate 78 multifractal approach, Zhou 2008). Instead of calculating the fluctuation function of the detrended variance using equation (3.4), a detrended covariance in each segment of size s was calculated as follows: 3.8 𝐹!,! 𝑠 = 1𝑠 [𝑌! 𝑖 − 𝑌!!!!!! 𝑖 ][𝑌! 𝑖 − 𝑌!! 𝑖 ] where Y(i) correspond to the integrated time series defined in equation (3.2), s refers to the segment length, and the superscripts o and f refer to oceanographic and fishery time series, respectively. Then, the qth order detrended covariance for non-zero q values was calculated as: 3.9 𝐹!,! 𝑞, 𝑠 = 1𝑁! 𝐹!! !!!!!!! Again, in an analogous way to equation (3.7) we expect that, in the presence of scaling, the following relation between two multifractal time series can be established: 3.10 𝐹!,! 𝑞, 𝑠 ∼ 𝑠!(!,!) It is clear that when the same time series is used in equation (3.9), which means superscript o equals superscript f, the above methodology reduces to the MF-DFA, i.e., the detrended covariance reduces to the detrended variance (Podobnik and Stanley 2008). The generalized Hurst coefficient defined in equation (3.7) for univariate multifractals can be expressed for a bivariate multifractal relationship as: 79 3.11 ℎ!,! 𝑞 = ℎ!,! 𝑞 + ℎ!,! 𝑞 /2 This latter equation means that the generalized Hurst exponent for two long-range dependent time series that are also long-range cross-correlated (ho,f), correspond to the average of the generalized Hurst exponents estimated individually for each multifractal oceanographic (ho,o) and fishery (hf,f) time series. 3.3 Results 3.3.1 Filtered sea surface temperature and wind stress time series The sea surface temperature time series shows a marked annual cycle, which is constant through the whole study period (1940-2013). This is clearly shown by the continuous wavelet global spectrum (Torrence and Compo 1998), where the annual cycle is marked as a continuous red band centered on a scale equal to 365 days (Fig. 3.5). Figure 3.5. Sea surface temperature wavelet power spectrum for 1940-2013 period. The Morlet wavelet was selected as the mother wavelet. The white dotted line indicates the cone of influence (COI), below which power spectrum coefficients should be interpreted with caution. A frequency band in dark red associated with the annual cycle is clearly visible. 80 Normalized sea surface temperature time series can be observed in Figure 3.6a. The reconstructed SST anomaly time series is shown in Figure 3.6b. Figure 3.6. a) Time series of normalized sea surface temperature (SST in °C) for the 1940-2013 period, and b) Wavelet reconstructed sea surface temperature representing variability between 128 and 4096 days (0.35 and 11.22 years, respectively). The period June 1997- June 1998 (highlighted with a orange bar) represents one of the two strongest El Niño events detected during the last century. Vertical lines separate periods of ten years. The filtered sea surface temperature time series (SSTF), used for the estimation of long-range dependence (Fig. A.1a, Appendix A), follows closely a Gaussian distribution (Figure A.1b and Figure A.1c), which was previously detected by Cummins and Masson (2014) for sea surface temperature (SST) anomalies in the same study area. The slope (-β) of the power spectrum plot (see Methods in Chapter 2) is equal to 0.63 which means that SSTF varies between -1<β<1, and as a consequence, it can be modeled as a fractional Gaussian noise time series (Seuront 2010). 81 On the other hand, normalized and reconstructed wind stress time series are shown in Figures 3.7a and 3.7b, respectively. The slope (-β) of the power spectrum plot for wind stress filtered time series (WSF) equals 0.72, which means that it behaves similar to a fractional Gaussian noise, however, a closer inspection to the quantile plot (Fig. A.2b, Appendix A), revealed that WSF shows clear deviations from Gaussianity. Figure 3.7. a) Time series of normalized daily mean wind stress (WS, m/s2) for the 1958-1998 period, and b) wavelet reconstructed wind stress representing variability between 128 and 4096 days (0.35 and 11.22 years, respectively). Vertical lines separate periods of ten years. 3.3.2 Long-range dependence of oceanographic and fishery time series According to the results obtained in the previous section, sea surface temperature time series can be modeled as fractional Gaussian noise (fGn time series) (in contrast to wind stress filtered time series) making possible the estimation of Hurst coefficients (H) using 82 the well-known maximum likelihood theory. Maximum likelihood estimation of H is possible only when the analyzed time series behaves as fGn (or summed fGn). Due to the specificity of this method, it cannot be used for a wide class of fractal time series. Accordingly it was not selected in the methodology section as the principal method for the estimation of LRD. In consequence, it was used here as an auxiliary methodology when fGn time series were obtained, enabling the estimation of Hurst coefficients that are independent of selected scaling ranges. Before the estimation of H using the DFA method, maximum likelihood estimates of Hurst coefficients (HML) were obtained for annual SSTF time series from 1940 to 2013 (Appendix A.3, Fig. 3.8). The intrinsic assumption that SSTF can be represented by a single LRD parameter that does not vary over time -which is valid only for long and exactly simulated LRD time series- cannot be sustained as expected (Fig. 3.8). However, relatively constant Hurst coefficients above a long-term mean were obtained for prolonged periods: 1940-1953, 1959-1974 and 1976-1986 (visualized as runs with similar values in Fig. 3.8). Those periods were interrupted during 1957 and 1975, years coinciding with a major El Niño event and a regime shift in the North-East Pacific, respectively. After 1996 and until 2013 the behaviour of H changed markedly from one year to the next; and again marked low H values were obtained for years dominated by anomalous oceanographic conditions: 1998 and 2011 (Fig. 3.8). For example, the lowest SST filtered daily values from the whole time series were detected during the major 1997-1998 El Niño event (Fig. A.4, Appendix A). Except for El Niño event of 1983, major oceanographic changes were detected in filtered sea surface temperature time series by marked decreases of Hurst coefficients calculated by the maximum likelihood method (HML) (Fig. 3.8). 83 Figure 3.8. Hurst coefficients estimated for filtered sea surface temperature (SSTF) time series during the 1940-2013 period. The horizontal dotted line indicates the mean of Hurst coefficients (0.9661). Vertical grey lines separate periods of ten years. Major El Niño events (1957, 1983, 1998) are demarked with filled red circles. Year 2011 (filled orange circle) coincides with an alternation of strong El Niño/La Niña events. Regime shifts detected in the North-East Pacific for 1988 (Hare and Mantua 2000) and 1975 (Overland et al. 2008) were highlighted with a filled blue circle. The average value of Hurst coefficients obtained for 74 years (1940-2013, 𝐻!"=0.9661) was used as a target or reference level to select the upper scaling range (smax) when H was estimated using Detrended Fluctuation Analysis (HDFA) (Fig. A.5, Appendix A). When a fractal time series behaves as fractional Gaussian noise (fGn), such as the filtered sea surface temperature time series, the Hurst coefficient equals the scaling exponent or LRD parameter α obtained from equation (3.5) (Taqqu et al. 1995, MacIntosh et al. 2013). This means that the upper scaling range (smax) was gradually shortened, starting from N/4 until a value equal to 𝐻!"=0.9661 was obtained. The selection of the scaling range using maximum likelihood estimates of H (HML) to obtain a representative value of the Hurst 84 coefficient using DFA (HDFA) is represented in the flow diagram on Figure A.5 (Appendix A). Figure 3.9 Root mean square fluctuation function F(s) for filtered sea surface temperature time series (SSTF) between 1940 and 2013. a) Scaling range (blue circles) extends between 16 and 1024 days (24-210, H=0.8760). A scaling break marked with an arrow is visible around 256 days. b) Scaling range (blue circles) extends between 16 and 80 days (24-26.32, H=0.9675). Lines were fitted to scaling ranges to estimate the Hurst coefficient. A scaling break was detected for sea surface temperature time series around 256 days (28; Fig. 3.9a). The Hurst coefficient obtained for the scaling range between 16 and 80 days estimated with Detrended Fluctuation Analysis method was equal to 0.9675 (Fig. 3.9b). Maximum likelihood estimates of H were not calculated for the filtered wind stress time series (WSF) because it showed deviations from a Gaussian distribution (Fig. A.2, Appendix 85 A). Instead, the Hurst coefficient was directly estimated using Detrended Fluctuation Analysis for the whole time series (N=14965). Two scaling ranges were detected (Fig. 3.10a,b), and the first, which extends between 16 and 120 days, was used for the estimation of the Hurst coefficient (H=0.7538). Figure 3.10 Root mean square fluctuation function F(s) for filtered wind stress time series (WSF) for the 1958-1998 period. a) First scaling range (grey line fitted to blue circles) extends between 16 and 120 days (24-26.90) (H=0.7538). b) Second scaling range (grey line fitted to red circles) extends between 512 and 3726 day (H=0.6563). For smooth pink shrimp total daily catch anomalies (TDCa) two clear scaling ranges were 86 detected. The first extended between 16 and 120 days (Fig. 3.11a) and the second between 121-642 days (Fig. 3.11b). Figure 3.11 Root mean square fluctuation function F(s) for smooth pink shrimp total daily catch (TDC) anomaly time series for the 1987-1998 period. a) First scaling range (grey line fitted to blue circles) extends between 16 days and 120 days (24- 26.90, H=0.8568). b) Second scaling range (grey line fitted to red circles) extends between 121 and 642 days (H=0.9681). The first scaling range coincides with that obtained when the Hurst coefficient was estimated for a much shorter time series of catches in Chapter 2 (H=0.8779, N=642); H obtained for this longer series is slightly lower (H=0.8568, N=2568). This result confirms 87 that long-range dependence is present in smooth pink shrimp total daily catches and it does not depend on the phase or period of the fishery. 3.3.3 Multifractality in oceanographic and fishery time series Generalized Hurst exponents (h(q), equation 3.7) for filtered sea surface temperature (SSTF) did not vary remarkably. Slowly decreasing exponents were obtained for statistical moments q=1 (h=0.9951), q=2 (h=0.9675), q=3 (h=0.9415), q=4 (h=0.9170). For this time series length (N=27010), bigger differences between generalized Hurst exponents are expected for multifractal time series (López and Contreras 2013). In fact, the difference between h(q)’s at q=1 and q=4 is less than 10%, and for this reason, SSTF time series should be considered as (mono)fractal. Similar generalized Hurst exponents were obtained for filtered wind stress time series (WSF, N=14965). Hurst exponents show minor differences for q=1 (h=0.7637), q=2 (h=0.7538), q=3 (h=0.7476), and q=4 (h=0.7431). The h(q) calculated for smooth pink shrimp total daily catch anomalies (TDCa) showed similar variability to that estimated for the oceanographic time series. Homogeneous exponents were obtained for all statistical moments: q=1 (h=0.8722), q=2 (h=0.8568), q=3 (h=0.8500) and q=4 (h=0.8443). As a consequence, TDCa is also classified as mono(fractal), and the estimation of multifractal exponents was not necessary (Fig. A.5b). 3.3.4 Bivariate long-range dependence between oceanographic and fishery time series The bivariate relationships between filtered sea surface temperature and total daily catch anomalies was tested considering the smallest scaling range detected for both time series independently (16-80 days). The estimated bivariate Hurst exponent (equation 3.10, Fig. 3.12) is ho,f = 0.8580 (N=2568). A close agreement was found between this latter value and 88 the theoretical value (ho,f = 0.8560) expected when bivariate long-range cross-correlation was tested between SSTF and TDCa using equation (3.11). This means that scaling ranges were well defined, and consequently the generalized Hurst exponent was accurately estimated (Fig. 3.12). The bivariate Hurst exponent ho,f between filtered wind stress time series and total daily catch anomalies equals 0.8315 (N=2568) (Fig. 3.13). Figure 3.12 Bivariate long-range cross-correlations between filtered sea surface temperature and smooth pink shrimp (TDC) anomaly time series for the 1987-1998 period. Scaling range (highlighted with filled circles) used to estimate the bivariate Hurst exponent extends between 16 days and 80 days (24- 26.32, ho,f=0.8580). 89 Figure 3.13 Bivariate long-range cross-correlations between filtered wind stress and smooth pink shrimp (TDC) anomaly time series for the 1987-1998 period. Scaling range (highlighted with filled circles) used to estimate the bivariate Hurst exponent extends between 16 days and 120 days (24- 26.90, ho,f=0.8315). 3.4 Discussion 3.4.1 Long-range dependence in oceanographic time series Long-range dependence with a marked Hurst coefficient (H=0.9675) was detected for sea surface temperature (SST) off the west coast off Vancouver Island (WCVI). This means that SST observations are correlated up to 80 days for the first scaling range (between 16 and 80 days) (Fig. 3.9b). Similar results were found for SST recorded at several lighthouse stations off the WCVI, including Amphitrite Point (Cummins and Masson 2014). These researchers detected decorrelation times that vary between 3 and 9 months. The first value (90 days) is 90 similar to the upper scaling limit (80 days) obtained here for the first scaling range of SSTF time series (3.9b). Sea surface temperatures from British Columbia lighthouses appear to co-vary with air temperature anomalies from the northern hemisphere over a range between 450 days (15 months) and 1825 days (5 years) (Freeland 1990). This suggests that the dynamics of sea surface temperature of the coastal zone is the result of an integration of processes acting on a basin to global scale, such as El Niño and/or Pacific Decadal Oscillations (Percival et al. 2008, Amos et al. 2015). In fact, significant co-variability between Amphitrite Point Lighthouse sea surface temperature and the Pacific Decadal Oscillation at low-frequency scales (periods >1.5 years) was reported by Cummins and Masson (2014). According to our results (Appendix A.3), the magnitudes of Hurst coefficients (H) for the SST series are affected by those large scale processes (Fig. 3.8), pulling down H towards a white noise process (for which H=0.5). However, the removal of major variability associated with low-frequency oscillations (section 3.2.3) buffered their distorting effect, and restricted the variability of H to a narrow band, as discussed in the section 3.3.2. For the wind stress (WS) time series, a clear scaling range (16-120 days) was detected. The H coefficient for this range is lower (H=0.7538) than that estimated for sea surface temperature. The strength of wind induced upwelling, during spring/summer shrimp fishing seasons, has increased from 1959 due to a rise in alongshore wind stress (Foreman et al. 2011). This increment could explain persistence detected in wind stress time series. An association between the coastal upwelling low frequency signal with the Pacific Decadal Oscillation (PDO) along the North Pacific coast was found by Macias et al. (2012). This reaffirms our hypothesis that long-range dependence detected for wind stress series could be explained, to 91 some extent, by large-scale atmospheric and oceanographic processes. For example, long-range dependence was estimated for the PDO index using different methods (including Detrended Fluctuation Analysis), and marked persistence in Hurst coefficients (H>0.5) was found at monthly and annual scales (Khaliq and Gachon 2010, Yuan et al. 2014). The long-range cross-correlation found between filtered sea surface temperature (SSTF) and TDC anomaly time series (ho,f=0.8580) has implications for the management of the smooth pink shrimp fishery. Under this perspective, long-range dependent (LRD) oceanographic time series are not only indications of large-scale atmospheric processes. It appears that LRD detected in SSTF time series is affecting the abundance and/or distribution of smooth pink shrimp populations off the west coast of Vancouver Island at seasonal scales (between 16 and 80 days, Fig. 3.12), or at least how the fishery responds to these variations. This result extends previous findings detected for this fishery, in which short-term effects of sea surface temperature (as water mass type) and wind stress were identified as oceanographic factors affecting shrimp availability at scales shorter than 15 days (Perry et al. 2000). 3.4.2 Variability and representativeness of Hurst coefficients Maximum likelihood Hurst coefficients (HML) obtained for filtered sea surface temperature (SSTF) during years 1957 and 1998 show low values in comparison to the time series mean (𝐻!"=0.9661), which coincides with strong El Niño events (Tully et al. 1960, Rasmusson and Carpenter 1982, Emery and Hamilton 1985, Hsieh et al. 1995, Wolter and Timlin 2011) (Fig. 3.8). In addition, marked changes in the trend of H occur in the middle of the 1970’s (1975) and end of the 1980’s (1988/1989) coinciding with regime shifts in the Northeast Pacific (Hare and 92 Mantua 2000, Overland et al. 2008) and strong El Niño events that started during the middle of 1975 and 1988 (Wolter and Timlin 2011). The regime shift of 1975 was detected using the same climatic index that detected the recognized regime shift of 1976/1977 (Hare and Mantua 2000). This index represents the mid-level upper atmospheric circulation and is known as Pacific North American pattern (PNA) (Overland et al. 2008). In addition, for 1975 a major change or regime shift in the North Pacific Ocean was detected which decrease the temperature down to 250 m after that year (Stephens et al. 2001). The lowest Hurst coefficient obtained for annual SSTF time series was for 2011 (H=0.835). This year marked the end of a rapid transition of strong El Niño/La Niña events (2010/2011) estimated using Amphitrite Point sea surface temperature (Hourston and Thomson 2013). These results show that although Hurst coefficients (HML) decreased during periods characterized by “anomalous” oceanographic conditions, they still fluctuate within a narrow band (0.83<HML<0.99) in comparison to the full range of possible H values (0<H<1). According to Fig. 3.8, Hurst coefficients always show long-range persistence (H>0.5, Taqqu et al. 1995), and its recovery -measured as the tendency to return to values closer or above historical mean (𝐻!"=0.9661)- is fast and occurs during the subsequent one or two years. The only exception for this general behaviour occurs after 1997-1998, years representing one of the two strongest El Niño events detected during the last century (Wolter and Timlin 2011). In fact, the lowest values of filtered sea surface temperature are concentrated within this latter El Niño period (Fig. A.4, Appendix A). Negative values of temperature were detected for the 1997-1998 El-Niño event because oscillations, which include El-Niño frequency, were previously removed for the elaboration of filtered sea surface temperature 93 time series. The original sea surface temperature clearly detected the increase in temperature for that period (Fig 3.6a). In summary, the average value of H maximum likelihood coefficients (𝐻!"=0.9661), used for selecting the scaling range when H was obtained by Detrended Fluctuation Analysis, represents the general long-range dependence behaviour of sea surface temperature, regardless of the increased variability shown after 1997-1998 (Fig. 3.8). 3.4.3 Selection of a long-range dependence over a short-range dependence model A null hypothesis behind this Chapter is that a long-range dependent (LRD) model is used to describe a first order short-range autoregressive (AR1) process. The model and process are switched, so we expect to reject the null hypothesis. It has been argued that autoregressive and long-range dependence models seem to fit equally well dependence between observations in time series of short/medium length which makes difficult the discrimination between a short-range and a long-range memory process (Maraun et al. 2004). Therefore, a quantitative approach is necessary for assessing the adequacy of these models. Traditionally, oceanographic and fisheries time series are modeled as short-range autoregressive processes (Hall and Manabe 1997, Pyper and Peterman 1998, Rudnick and Davies 2003); however, as monitoring continues and time series length increases a growing amount of evidence against this hypothesis is accumulated. The probabilities of erroneously using a long-range dependent (LRD) model to describe a short-range autoregressive (AR1) process were calculated using Monte Carlo simulation (Appendix A, Fig. A.6). As observed on Figure A.6 (Appendix A), this probability depends on time series length. At lengths equal to shrimp total daily catch, wind stress and sea surface temperature time 94 series, the probabilities of rejecting the null hypothesis are high and equal to 0.8977, 0.9905, 0.9999, respectively (Fig. A.6, Appendix A). In the worst case scenario, when using the shortest time series which corresponds to pink shrimp total daily catch (N=2568), the probability of falsely describing an autoregressive process using a long-range dependent model (using high LRD and AR1 coefficients, Fig. A.6, Appendix A) equal to 10%. This means that the choice of using a LRD model to fit smooth pink shrimp total daily catch time series was adequate and supported by a 90% of chance of correctly selecting the model. It is important to remark that this result does not mean that when analyzing pink shrimp total daily catch time series to estimate the Hurst coefficient, there is a 10% probability that long-range dependence is not present (falsely detecting long-range dependence). It only means that 10% of the time we used the wrong model to describe a process, and longer time series are required to confirm with a probability higher than 90% the correct selection of the model (according to the Spectral Density Test, Fig. A.6, Appendix A). It is clearly observed on this latter figure that when pronounced long-range dependent (H) and autoregressive (AR1) coefficients are used in our simulations (both are considerably greater than 0.5), the time series length needed to reduce the wrong selection probability from 10% to a value close to 1% gets close to N=15,000 observations. These results mean that we are confident in using a LRD model to estimate Hurst coefficients for shrimp total daily catches time series, and specially for sea surface temperature and wind stress, because for both of them the probabilities of using the wrong model fluctuate around 1% or less (Fig. A.6, Appendix A). Future studies should also test the ability of alternative long-range dependent models to fit empirical observations as 95 recommended by Seuront et al. (2015). Long-range dependence has been detected in oceanographic time series from different marine ecosystems worldwide. Smooth pink shrimp fishery management areas are located in a complex and dynamic oceanographic transition zone, strongly influenced by superposed climate cycles. Detrended fluctuation analysis (DFA) was used to detect and quantify long-range dependence in sea surface temperature; wind stress and total daily catch (TDC) time series from the west coast off Vancouver Island. Oceanographic and fishery time series showed no signs of multifractal patterns of variability, and accordingly were classified as monofractal. Hurst coefficients estimated for sea surface temperature were affected by El Niño events and regime shifts; the lowest values of H were obtained for years characterized by “anomalous” oceanographic conditions. The same scaling range detected in Chapter 2 for total daily catches (between 16 and 120 fishing days) was detected in the present chapter when a time series of catches four times longer was analyzed. This result confirms that long-range dependence (LRD) is present in total daily catches, and that the estimation of LRD was possible despite the presence of missing fishing days. Bivariate long-range dependence was detected among TDC and sea surface temperature (scaling range between 16 and 80 fishing days), and among TDC and wind stress (scaling range between 16 and 120 fishing days). These results may imply that sea surface temperature and wind stress are affecting the abundance and/or distribution of smooth pink shrimp populations at intraseasonal scales, or at least how the fishery responds to these variations. On the other hand, these results indicate that having quantified the variability of total daily catches (TDC) at small temporal scales (16 fishing days), it is 96 possible to infer at any time during the progress of the fishing season the variability of smooth pink shrimp catches up to 120 fishing days. However, analyses focused on total daily catches (TDC) conducted in Chapter 2 using the structure function (SF) methodology, and those conducted in Chapter 3 using the Detrended Fluctuation Analysis (DFA) are not conclusive. Results obtained in Chapter 2 suggest that total daily catches are multifractal, however, when a long time series of catches was analyzed in Chapter 3, weak evidence of multifractality was detected. In this context, a new hypothesis was tested in Chapter 4 when early warning indicators were constructed. Under this new hypothesis framework multifractality can vary over time, which could explain differences in results obtained in Chapter 2 and Chapter 3. For that reason, in Chapter 4 indicators were constructed quantifying the multifractal dynamics of total daily catches within rolling windows which moves across the history of the fishery. 97 Chapter 4: A self-similar early warning indicator accounting for long-range dependence and intermittent variability in Pandalus jordani fishery catches from the west coast of Vancouver Island 4.1 Introduction Two main factors have been proposed, and to a great extent demonstrated, to be responsible for dramatic changes of the status of fish populations: high exploitation rates (Myers et al. 1997), and the combination of overfishing and adverse environmental conditions (MacCall 2011). However, it is clear that anthropogenic pressure (overfishing) plays a major role in moving a fishery towards an undesirable state (depleted or collapsed, Hammer et al. 2010). The transition between a desirable and undesirable state occurs after surpassing some threshold after which the recovery of the fishery can be extremely long and difficult (Hutchings 2000, Hutchings & Reynolds 2004) and in some cases even not possible after a major reduction in exploitation rates (Worm et al. 2009, Petitgas et al. 2010). The economic costs of restoring overfished stocks are by far greater than adopting proactive management measures in order to avoid a fisheries collapse (Gaines and Costello 2013). On the other hand, the recovery of a fishery depends on a combination of factors, not all of them controllable by fisheries managers (Wakeford et al. 2009), and for that reason, avoiding stock degradation become a highly desirable objective. An increased amount of research has been conducted on analyzing temporal changes in exploited and non-exploited ecosystems, focusing on detecting large scale abrupt changes (regime shifts) and critical thresholds over time, in order to give managers the chance to intervene and avoid the collapse of populations (Batt et al. 2013, Carpenter et al. 2014). 98 The detection of these thresholds has been conducted using different classes of early warning empirical indicators proven to be useful in detecting major changes in complex systems, including aquatic ecosystems (Dakos et al. 2008, Biggs et al. 2009). Indicators used to detect in advance ecosystem changes (leading indicators) can be classified into three main categories according to the statistical time series properties they are trying to detect (Dakos et al. 2012). The first category of indicators is based on the detection of an increase in time series memory before a major ecosystem change takes place by analyzing its autocorrelation structure in time and frequency domains (Kleinen et al. 2003, Dakos et al. 2008, Dakos et al. 2012). The second category uses some estimator of variability (e.g. variance, skewness) in an attempt to detect their increase through time, as a major change in the ecosystem structure becomes imminent (Carpenter and Brock 2006, Guttal and Jayaprakash 2008, Seekell et al. 2011). The third category focuses on the detection of time series short jumps (called flickering) into an alternative ecosystem state before the complete shift or critical transition takes place (Dakos et al. 2013). Each class of indicators attempts to detect increased ecosystem sensitivity to perturbations, which appears to be characteristic close to the point where a major change takes place. It appears that only a few cases of their application have been reported in the fisheries literature. For example, spatial variance indicators have been used to detect collapses of crustacean fisheries (Litzow et al. 2013), and time series variance indicators have been applied to detect major shifts in simulated time series obtained from fisheries food-web models (Biggs et al. 2009). In general terms, proposed empirical indicators only represent limited aspects of time series dynamics, such as short-range memory structure (e.g. autocorrelation) or some estimator of variability (e.g. variance, kurtosis, skewness). 99 In the present study, taking into account changes in the temporal dynamics of indicators exhibited previous to major ecosystem reorganizations (Scheffer et al. 2012, Dakos et al. 2013), we constructed fisheries time series indicators based on the fractal pattern of variability detected for the smooth pink shrimp (Pandalus jordani) fishery from the west coast of Vancouver Island (Chapters 2 and 3; Montes et al. 2012). This pattern implies the presence of a slow decaying autocorrelation function (long-range dependence, LRD) and high variability (intermittency) in time series of catches. The superposition of LRD and non-Gaussianity phenomena (due to the presence of rare or intermittent observations) in ecosystems can be expressed as self-similarity (H, Franzke et al. 2012), a basic property of fractals (introduced in Chapter 1). Under this framework, self-similarity integrates time series long-range dependence (estimated by the LRD parameter d) and variability or intermittency (estimated by the stability exponent α) into a single parameter called the Hurst coefficient H (Taqqu and Teverovsky 1998, Frankze et al. 2012). LRD is synonymous with persistence, and intermittency (non-Gaussian jumps) refers to large changes that can occur abruptly, forcing a time series dominated by this latter effect to exhibit sharp discontinuities. LRD and intermittency are known as Joseph and Noah effects, respectively (Mandelbrot and Wallis 1968), and arise simultaneously pulling ecosystem dynamics to move into opposite directions (Frankze et al. 2012). The performance of LRD estimators under the assumption of Gaussianity is well studied (Taqqu et al. 1995, Taqqu and Teverovsky 1998, Gençay et al. 2002, Robbertse and Lombard 2012), however, under a non-Gaussian framework (non-Gaussian jumps) the 100 performance of the LRD estimator (d) and intermittency estimator (α) varies greatly, and caution should be taken in the presence of trends and noise (Franzke et al. 2012). In this Chapter, the first goal was to construct early warning indicators based on the evidence of multifractality detected in smooth pink shrimp TDC time series in Chapter 2. In contrast to analyses conducted in Chapter 2 and Chapter 3, in this chapter multifractality was quantified within rolling windows to test the hypothesis that multifractal dynamics can vary through time. The second goal was to verify the usefulness of proposed indicators to track in-season dynamics of smooth pink shrimp TDC and detect in advance temporal changes underlying this fishery, using historical biomass as an indicator of the status of these smooth pink shrimp populations. The third goal was to compare the performance of self-similar indicators with classical statistical time series indicators proposed in the literature to detect in advance changes in aquatic ecosystems under exploitation. 4.2 Methods First, a well-known long-range dependence model (FARIMA model, Taqqu and Teverovsky 1998) was used to simulate LRD time series with different intermittency levels (α values). The maximum number of TDC missing observations detected for a single fishing season of the smooth pink shrimp fishery was introduced into the simulated time series and then Hurst coefficients for gappy data, estimated using two different methods, were compared with theoretical H values. In this context, uncertainties and bias in H due the presence of missing observations were estimated using simulated data, and then using 95% confidence intervals the best method was selected and applied for the estimation of Hurst coefficients of smooth pink shrimp total daily catches localized in time (known as local Hölder exponents, LHE’s). 101 Before the estimation of LHE’s, marked changes in fishing effort from one fishing season to the next were detected and removed using continuous and discrete wavelet transforms (Percival and Walden 2000). 4.2.1 Estimation of Hurst coefficient confidence intervals for simulated gappy data. Total daily smooth pink shrimp catch time series for the 1987-1998 period (N=2568, 12 years), which contains 16.5% of missing observations, was analyzed. Missing observations for each fishing season, constructed in the same way as in Chapter 2, were calculated. They fluctuate between 1% (for 1995 and 1996 fishing seasons) and 32% (1990 fishing season) and consequently no attempt to interpolate these values was made. Missing values were replaced by zeros, which tends to underestimate LRD coefficients only at short time lags. However, this effect is well studied for time series of the order of 10,000 data points (Wilson et al. 2003). In this context, we preferred to quantify the effect that replacing missing values with zeros into pink shrimp TDC time series can have on LRD coefficients using simulated time series. We simulated a LRD time series with intermittent variability using a FARIMA model with symmetric alpha stable (SαS) innovations (Stoev and Taqqu 2004, 2005). FARIMA (p,d,q) models extend short-range ARIMA (p,d,q) models replacing the differencing integer operator d used in ARIMA models with a fractional real number (Beran et al. 2013). The same operator d has different meanings in these two models. In short-range ARIMA models d represents the degree of differentiation needed to remove linear or seasonal trends in order to achieve a stationary condition (variance and mean approximately constant over time). In FARIMA models d represents the long-range memory parameter which fluctuates 102 between 0<d<1-1/α, and α is the degree of intermittency which fluctuates between 1<α<2 (Stoev and Taqqu 2004, Beran et al. 2013). In summary, a FARIMA model with symmetric alpha stable (SαS) innovations incorporates four components into a single model: i) a short-range autoregressive (AR) component represented by parameter p, ii) a long-range fractionally differenced or I component (represented by the non-integer parameter d), iii) a correlated random component (also called moving average MA, represented by parameter q) and iv) a highly intermittent component (alpha stable, represented by parameter α). In the simulations, the autoregressive coefficient was fixed to a value commonly found in fisheries time series (p=0.65, Pyper and Peterman 1998) and the moving average coefficient (q) was set to zero. The length of the simulated time series was equal to the length of the rolling window (rw) used to estimate self-similar parameters (rw=1284, section 4.2.3). Time series replicates were obtained for different H=d+1/α values, starting from d=0.1071 and α=1.6 to d=0.5 and α=2. For each value of α (1.6, 1.7, 1.8, 1.9, 2.0), six different LRD parameters d were estimated. To account for the effect of missing observations, we adopted a very conservative approach in which the maximum percentage of missing values detected for TDC time series using a rolling window equal to 1284 was considered. Missing values varied between 7.5 and 25%, and this latter percentage of missing observations was introduced into the simulated time series. We also introduced a sequence of 17 continuous zeros into the simulated time series, which was the largest segment of continuous missing observations found for the smooth pink shrimp TDC time series. 103 Finally, replicates (R) were obtained for 1,000 resamples using the maximum entropy bootstrap (ME bootstrap, Vinod 2006, Vinod and López-de-Lacalle 2010, Barbosa 2011) which preserves the LRD structure of the original time series. Using this approach, confidence intervals (CI) for d and α estimated at a 95% level were obtained. These CI’s reflect the uncertainty of self-similar parameters estimated for a time series equal to 1284 observations considering the highest fraction of missing data (25%) that can be detected in each rolling window. 4.2.2 Complexity of TDC time series and removal of unwanted oscillations In general terms smooth pink shrimp TDC exhibited complex dynamics, characterized by several embedded process acting at different time scales, which is a common pattern for fisheries time series using the continuous wavelet transform (Menard et al. 2007, Rouyer et al. 2008). Wavelets are functions capable of narrowing to focus on high frequency or widening to focus on low-frequency components of a time series. They are parameterized by a translation parameter b (-∞<b<∞), which denotes the location of the time series being analyzed, and the scale parameter a (a>0) which corresponds to the magnification (scale or resolution) at location b (Addison 2002, Mallat 2009). A wavelet base can be constructed from a unique function ψ(x) according to: 4.1 𝜓!,! 𝑥 = 𝜓 𝑥 − 𝑏𝑎 where ψ is a single function known as mother wavelet and x represents the time series axis. Under this framework, a continuous wavelet transform (CWT; Torrence and Compo 1998, Percival and Walden 2000) is expressed as: 104 (4.2) 𝑊[ℎ] 𝑎, 𝑏 = !! 𝜓!,!∗!!! 𝑥 ℎ(𝑥)𝑑𝑥 where ψ* denotes the complex conjugate of ψ x . A continuous wavelet transform can be thought as a cross-correlation between a time series h(x) and a set of wavelets ψ*((x-b)/a). These latter functions are dilated (stretched) and translated (shifted) onto the time series allowing the extraction of time-dependent amplitude changes. Wavelet coefficients W[h](a,b) represent the contribution of the different scales a at different time positions b. This means that coefficients are obtained from correlating segments of a time series with translated and dilated versions of the selected wavelet function. The term 1 a is a normalization factor that allows a comparison for different scales (Mallat 2009). When these equations are applied to the smooth pink shrimp TDC time series, patterns of variability at multiple time scales are apparent (Fig. 4.1). Those scales demarcated with black contours encompass regions significantly different from an autoregressive (AR) process with lag-1 coefficient equal to 0.65 (Fig. 4.1). This means that a more complex model, different from the AR1 model used to quantify the correlation structure of catches distant by one fishing day, is necessary to describe the temporal structure of daily catches. Red zones indicate a similarity between TDC time series and the mother wavelet (Morlet), which is frequently used to detect periodic patterns of variability and is well localized in scales (Cazelles et al. 2008). Definitions of variables and parameters used in the present chapter are summarized in Table 4.1. 105 Table 4.1. Definitions of variables and parameters used in Chapter 4. Symbol Definitions LRD Long-range dependence H Hurst coefficient d Fractional difference parameter, long-range dependence parameter 1/d Inverse long-range dependence parameter α Stability exponent, degree of intermittency or Lévy index AR, MA Autoregressive and moving average components or processes ARIMA Autoregressive integrated moving average model FARIMA Fractionally differenced autoregressive moving average model p,q Autoregressive and moving average parameters I Integer difference parameter for ARIMA models LHE Local Hölder exponent 1/LHE Inverse local Hölder exponent SαS Symmetric alpha stable innovations TDC Total daily catches 106 Symbol Definitions R Number of replicates or resamples ME Maximum entropy bootstrap CI Confidence intervals CWT Continuous wavelet transform ψ Mother wavelet ψ∗ Complex conjugate of mother wavelet a,b Scale and translation parameters from defined wavelet function h(x) Time series axis on the wavelet plane W[h](a,b) Wavelet coefficients at scale a and position b 1/ a Normalization factor COI Cone of influence TDC Total daily catch TDE Total daily fishing effort PFMA Pacific Fisheries Management Areas WCVI West coast off Vancouver Island 107 Symbol Definitions ML Maximum likelihood QQ Quantile method KG Kettani and Gubner estimator or method ρn(1) Sample autocorrelation at lag 1 rw Rolling windows AR(1) Autoregressive coefficient at lag 1 SR Spectral ratio sd Standard deviation LRP Limit reference point SSB Spawning stock biomass CSD Critical slowing down 108 Figure 4.1. Wavelet power spectrum for smooth pink shrimp total daily catch (TDC) time series for 1987-1998 fishing seasons. Vertical axis represents scales in days. Shaded area between 128 and 512 days contains most of the total daily catch time series variance. The Morlet wavelet was selected as the mother wavelet. Black contours encompass regions for which the lag-1 autoregressive coefficient (AR) is significantly different from 0.65. White dotted lines indicates the cone of influence (COI), below which the power spectrum coefficient should be interpreted cautiously. Periodic patterns of variability and abrupt changes can affect the estimation of long-range memory exponents and therefore should be removed (Shen et al. 2007). Periodicities are visible at small scales (concentrated between 2 and 8 days, Fig. 4.1), possibly associated with the mean duration of fishing trips (e.g. Chapter 2). After fishing day 2000, which correspond to 1995 fishing season, this periodicity disappears. It could be related to changes in fishing trips dynamics, because starting from this year 60% of smooth pink shrimp catches comes only from an area around Barkley Sound which is closer to the coast (Perry et al. 2000), and belongs to Pacific Fisheries Management Area 123 (Fig. 2.1). 109 Anyways, those patterns were not removed because they only affect time series variability at short time lags. On the other hand, at large scales, most of the TDC time series variance is heterogeneously distributed between 128 and 512 days (shaded area on Fig. 4.1), inducing non-stationarity, which means changes in the time series mean over time. The non-stationary dynamics of the TDC time series can be better observed in Figures 4.2a,b. In Figure 4.2b, breaks are visible at 1988, 1991, 1994, 1995 and 1996 fishing seasons when catch changes markedly. Those changes in catches appear to be associated with marked changes in fishing effort for those years, as reported by Martell et al. (2000), and also visible in Figure 4.2b. Deterministic breaks, also called step-like changes, affect the time series mean and can bias the estimation of Hurst exponents toward larger values (Shen et al. 2007, Rust et al. 2008, Rea et al. 2011). In the smooth pink shrimp TDC time series, those breaks were detected using discrete wavelet transform (Stoev et al. 2005, Shen et al. 2007), sampling the minimum number of CWT coefficients to represent time series variability at the 128-512 scale range (Gencay et al. 2002). As our ultimate goal was to detect early changes underlying the smooth pink shrimp fishery, the Haar wavelet basis was selected as the mother wavelet because this basis is more localized in time than on scales in comparison to Morlet wavelet (Cazelles et al. 2008, Lovejoy and Schertzer 2012). Detrending the TDC time series using wavelets coefficients is equivalent to using a band-pass filter (Gencay et al. 2002) to remove unwanted oscillations in fisheries catch data (e.g. Ménard et al. 2007). 110 Figure 4.2. a) Time series of smooth pink shrimp total daily catch (TDC) for 1987-1998 fishing seasons, b) catch level shifts from 128 to 512 days (step-like changes) were calculated using Haar wavelets. Total daily fishing effort (TDE in trawl hours, green line) follows major step-like changes in TDC, and c) TDC residuals obtained after subtracting step like changes in catches obtained in (b). Vertical grey lines separate annual fishing seasons. After detrending the smooth pink shrimp TDC time series, marked changes in catches associated with fishing effort dynamics fluctuating around 1 calendar year (between 128-512 days) were removed. In addition, the artificial periodicity associated with the stitching of fishing seasons (214 fishing days, Chapter 2) disappears. The residual time series, free of these unwanted oscillations, can be observed in Figure 4.2c. TDC residuals were used in the estimation of early warning signals. 111 Non-Gaussian fluctuations in TDC time series make the estimation of Hurst coefficients not possible using well known methods such as dispersional analysis (Caccia et al. 1997) or maximum likelihood methods (Deriche and Tewfik 1993, Pilgram and Kaplan 1998, Robbertse & Lombard 2012) for which Gaussianity is an assumption. Extreme small and large observations are frequently observed in the TDC time series causing intermittency, which is reflected as empirical distributions with heavier tails than a Gaussian distribution (Fig. 4.3). Intermittency and long-range dependence were parameterized to construct a single early warning index as explained in the following section. Figure 4.3. Histograms for the smooth pink shrimp total daily catch (TDC) fishing seasons from 1987 to 1998. 112 4.2.3 Estimation of Hurst coefficient, LRD and intermittency parameters The Hurst coefficient (H), which integrates long-range dependence and time series intermittency into a single measure, was used to detect early changes of TDC time series. H was calculated according to equation 4.3 (Weron et al. 2005, Franzke et al. 2012) as: 4.3 𝐻 = 𝑑 + 1/𝛼 where d corresponds to the long-range dependence (LRD) parameter, and α to the intermittency parameter. H was estimated using two different approaches. First, LRD (d) and intermittency (α) parameters were calculated independently, and added according to equation (4.3) to obtain estimates of the Hurst exponents. The long-range dependence parameter d, (Taqqu and Teverovsky 1998) was calculated using maximum likelihood (ML, Haslett and Raftery 1989) using the fracdiff R package (Fraley et al. 2006) from the R software (R Development Core Team 2005, version 3.0.2). The intermittency parameter α was calculated using the portes R package (Mahdi et al. 2014) according to the quantile method (QQ method) described in McCulloch (1986). This method, which was used to calculate H after the estimation of d and α using the R packages, is called from here as the QQ-ML method. A second method to estimate Hurst coefficients was used, which does not require the estimation of d and α separately (eq. 4.3), is easy to calculate (eq. 4.4) and provides estimates much faster than using the QQ-ML method. This second method, known as the KG estimator (Kettani and Gubner 2006), which uses only the autocorrelation coefficient estimated at the first time series lag to estimate H, was calculated as: 113 4.4 𝐻 = !! [ 1+ log (1+ 𝜌!(1) ] where ρn(1) denotes the sample autocorrelation at lag 1. This estimator, initially developed only for Gaussian processes (e.g. fractional ARIMA time series, Kettani and Gubner 2006), was recently proven to be robust against different types of noise. This means that the KG estimator can also be used to estimate the Hurst coefficient in time series with a heavy tail distribution, i.e. for a time series with some degree of intermittency (Sheng et al. 2011). Moreover, the KG estimator shows the best tracking performance among twelve H estimators when applied to processes characterized by time-varying Hurst estimators (Sheng et al. 2012a). These latter processes are known as multifractional processes where the single and constant Hurst exponent is replaced by a family of different exponents which estimates how self-similarity varies through time (Sheng et al. 2012a,b). It is reasonable to expect that the Hurst exponent for the smooth pink shrimp TDC time series will vary in response to anthropogenic, environmental, or other drivers. As a consequence, H cannot be assumed to be constant. In addition, it could vary through time in response to, for example, fisheries management actions or changes in population biomass. We adopted the approach used by Dakos et al. (2008, 2012) to detect early warning signals on a localized manner using rolling windows (rw) and linear trends. We used rolling windows (rw) of half the size of the TDC time series (rw=1284) to estimate Hurst exponents localized in time, which are called from here as local Hölder exponents (LHE). Using this approach, the first local Hölder exponent, calculated using the TDC observations from 1987 to 1992, is compared with the TDC obtained during the first day (April 1) of the 1993 fishing season and so on until the end of the TDC time series. 114 Changes in linear trends in the dynamics of indicators extracted from equation (4.3) (H, d and α) were used to detect and represent early warnings of major changes underlying this fishery. Local slopes were calculated to detect changes in rates of change of the indicators, and as a consequence, of linear trends. Local slopes fluctuating around zero can be associated with indicators following a linear trend (Takalo 1995, Maraun et al. 2004). On the contrary, slopes which deviate from zero denote breaks of linear trends and changes in the dynamics of indicators. Local Hölder exponents (LHE’s) were obtained using equation (4.3) and equation (4.4) and compared according to bias and uncertainty obtained from simulated time series. A trade-off between bias -measured as the difference between the average value of exponents and theoretical values-, and uncertainty -measured as the range of values within which theoretical values are contained-, was evaluated to select the most appropriate estimator. 4.2.4 Comparison between LHE’s and classical statistical indicators LHE’s proposed in this chapter were compared with statistical indicators for their ability to detect in advance major changes underlying the smooth pink shrimp fishery. Many statistical indicators have been proposed for detecting ecological changes (Guttal and Jayaprakash 2008, Dakos et al. 2011, Scheffer et al. 2012) but only those which have used fish catches and data obtained from fisheries models to detect early changes in aquatic ecosystems (Biggs et al. 2009, Litzow et al. 2013) were considered in this study. In this context, four different statistical indicators were selected. The first two indicators quantify time series memory properties: AR(1) coefficient as a measure of short-range memory, and spectral ratio (SR), which measures the ratio between spectral density at low 115 frequency (0.05) to the spectral density at high frequency (0.5) (Dakos et al. 2012). Then, in order to detect the rising variability usually found in fish catches before ecosystem collapses (Hsieh et al. 2006), skewness and kurtosis, defined as the standardized third and fourth statistical moments respectively, were used (Biggs et al. 2009, Carpenter et al. 2011, Litzow et al. 2013). When classical statistical indicators were calculated, rolling windows (rw) equal to half the length of the TDC time series (rw=1284) were used to match the same size of windows for the estimation of local Hölder exponents. For classical statistical indicators, missing TDC values were linearly interpolated, a logarithmic transformation (log TDC+1) and standardization (𝑇𝐷𝐶 −𝑚𝑒𝑎𝑛(𝑇𝐷𝐶)/𝑠𝑑(𝑇𝐷𝐶)) was conducted on the TDC time series (Fig. 4.4), as recommended by Dakos et al. (2012) for time series with values close to zero and extreme observations, respectively. Due to the extreme variability exhibited by the TDC series, the standardized TDC time series only approximates a Gaussian distribution (Fig. 4.4). 116 Figure 4.4. a) Standardized smooth pink shrimp total daily catch (TDC) time series from 1987 to 1998. Vertical grey lines separate annual fishing seasons. b) Histogram of standardized TDC time series. 4.3 Results 4.3.1 Hurst coefficient confidence intervals for simulated gappy data The comparison between QQ-ML and KG estimators obtained for simulated gappy time series shows that, in general, the KG estimator performs better than the QQ-ML estimator in terms of bias. In almost all cases, the KG estimator has an average value closer to the theoretical H value (horizontal red line) than does the QQ-ML estimator (Fig. 4.5 to Fig. 4.9). 117 Figure 4.5. Uncertainty and bias in the estimation of H from simulated gappy time series using the quantile/maximum likelihood method (QQ-ML) and the KG method as described in the text for α=1.6 and six increasing long-range dependence levels d. Red horizontal lines mark theoretical H value. Experimental H values observed in boxplots were calculated from replicates obtained using ME bootstrap. For low to medium values of LRD and all intermittency levels (top row for each α plot) which means from d=0.1071 for α=1.6 to d=0.2857 for α=2.0 (equivalent to H=0.7321 to 0.7857, respectively), the average KG estimator occurs above the theoretical reference line. 118 Figure 4.6. Uncertainty and bias in the estimation of H from simulated gappy time series using the quantile/maximum likelihood method (QQ-ML) and the KG method as described in the text for α=1.7 and six increasing long-range dependence levels d. Red horizontal lines mark theoretical H value. Experimental H values observed in boxplots were calculated from replicates obtained using ME bootstrap. In contrast, for intermediate to high values of LRD and all intermittency levels (bottom rows for each α plot), which means from d=0.2679 for α=1.6 to d=0.5 for α=2.0 (equivalent to H=0.8929 to H=1.000, respectively) (Fig. 4.5 to Fig. 4.9), the KG estimator is at or below the theoretical reference line. This means that, when applied to the smooth pink shrimp TDC time series, the KG estimator will slightly overestimate low and underestimate high values of LRD, respectively, in comparison to the results obtained using the QQ-ML method. 119 Figure 4.7. Uncertainty and bias in the estimation of H from simulated gappy time series using the quantile/maximum likelihood method (QQ-ML) and the KG method as described in the text for α=1.8 and six increasing long-range dependence levels d. Red horizontal lines mark theoretical H value. Experimental H values observed in boxplots were calculated from replicates obtained using ME bootstrap. When comparing the effect of low and high levels of intermittency in the estimation of intermediate to high values of H using the KG method, we detected that under extreme conditions, which means no intermittency (Gaussianity, α=2) or high intermittency (α=1.6) (Taqqu and Teverovsky 1998), similar results were obtained (Figures 4.5 and 4.9), which supports the conclusion that the KG estimator will not be affected by extremely low or high variability in the TDC time series. 120 Figure 4.8. Uncertainty and bias in the estimation of H from simulated gappy time series using the quantile/maximum likelihood method (QQ-ML) and the KG method as described in the text for α=1.9 and six increasing long-range dependence levels d. Red horizontal lines mark theoretical H value. Experimental H values observed in boxplots were calculated from replicates obtained using ME bootstrap. Both estimators, QQ-ML and KG, performed badly in terms of uncertainty, because the range of H values obtained using both estimators contain the theoretical value in only a few cases (Fig. 4.5 to Fig. 4.9). Note that the KG estimator is less variable than the QQ-ML estimator, which means that although the calculated values do not match the theoretical H coefficients, the individual and average values are most of the times closer to the theoretical values when the KG estimator was used. In few cases, for example when the Levy index α equals 1.9 and d 121 is greater than 0.4, the QQ-ML method performs better than the KG method (Fig. 4.8). This unbiased property of the KG method is useful for detecting changes in the TDC time series, because it can detect the combined effect that individual changes in LRD and intermittency have on the dynamics of the smooth pink shrimp TDC time series. Figure 4.9. Uncertainty and bias in the estimation of H from simulated gappy time series using the quantile/maximum likelihood method (QQ-ML) and the KG method as described in the text for α=2.0 and six increasing long-range dependence levels d. Red horizontal lines mark theoretical H value. Experimental H values observed in boxplots were calculated from replicates obtained using ME bootstrap. The QQ-ML method provides information about intermittency (α) and long-range dependence (d) independently (equation 4.3). However, the estimation of α is known to be difficult (Taqqu 122 and Teverovsky 1998), and possibly was the main contributor in obtaining biased and uncertain H values using this method (Fig. 4.5 to Fig. 4.9). Consequently, the intermittency parameter values (α) should be interpreted with caution. Based on this result, it can be concluded that the effect of missing values on the TDC time series can contribute to the uncertainty in H observed for different intermittency levels, however, they are not an obstacle to detecting changes in time series properties. As the KG estimator outperforms the QQ-ML estimator in quantifying self-similarity when time series are governed by LRD and intermittent dynamics, local Hölder exponents estimated by the KG method (equation 4.4) were considered more appropriate for use as an early warning indicator. 4.3.2 Local Hölder exponents, LRD and intermittency parameters 4.3.2.1 Early detection of changes in total daily catches (TDC) The fluctuations of local Hölder exponents (LHE) and local d (LRD coefficient) in relation to changes in TDC residuals are shown in Figures 4.10 and 4.11, respectively. In general, a similar pattern of variability is observed for LHE’s and local d, characterized by a continuous increase in magnitude, which means increasing self-similarity (Fig. 4.10) and increasing long-range dependence (Fig. 4.11) across the TDC residual time series. 123 Figure 4.10. Variability of local Hölder exponents (LHE’s, increasing red line) and total daily catch residual time series (black lines) for 1993-1998 fishing seasons. The most remarkable issue that can be observed is the marked increase in the rate of change of the local d (slope) at the beginning of 1995 fishing season. Then, before the end of the 1995 fishing season, the rate of increase of local d’s changed back to levels similar to those observed before the 1995 fishing season (Fig. 4.11). Changes in the rates of increase of local d’s were confirmed with a local slope plot (Fig. 4.12, Takalo 1995, Maraun et al. 2004, Chapter 2). In this plot, local slopes which deviate substantially from zero denote changes in the growth rates of local slopes. 124 Figure 4.11. Variability of local d’s (long-range dependence parameter, increasing red line) and total daily catch residual time series (black lines) for 1993-1998 fishing seasons. Figure 4.12. Local slopes estimated for 1993-1998 fishing seasons using the method described by Takalo (1995) and Maraun et al. (2004). A window equal to 91 data points was used for the estimation of d’s (long-range dependence). Major changes in slopes are visible for the 1995 fishing season. 125 The fluctuation of the local α (intermittency parameter or Lévy index) shows an increasing linear trend, which extends from 1993 up to early in the 1994-fishing season (Fig. 4.13). After that date and until the middle of the 1995 fishing season the intermittency parameter shows a zero slope because α surpass the threshold value of 2. An intermittency parameter of 2 means that the analyzed time series can be described by a log-normal distribution (Seuront 2010). As α decreases, the high values of the TDC time series do not dominate as much as for larger values of α. These dynamics can be observed for 1996, 1997 and 1998 fishing seasons (Fig. 4.13), characterized by relatively minor deviations from the mean TDC. Figure 4.13. Variability of local α (intermittency parameter, red line) and total daily catch residual time series (black line) for 1993-1998 fishing seasons. 126 4.3.2.2 Early detection of changes in vulnerable biomass Historical fluctuations of self-similar indicators estimated for TDC have been compared with changes of vulnerable biomass across the history of this smooth pink shrimp fishery. Vulnerable biomass, defined as the proportion of biomass accessible to the fishing gear at the beginning of each fishing season, was obtained from Martell et al. (2000). We are aware that vulnerable biomass is greater than biomass indices obtained from shrimp surveys starting from the 1987 fishing season, however, both biomass indicators follow approximately the same trends after this latter year (Martell et al. 2000). We included fluctuations of total daily catches on the same graphs to present an integrated view of changes in this smooth pink shrimp fishery. Inverse local Hölder exponents (1/LHE) estimated using the KG method are presented for a better visualization of early changes of exponents in relation to biomass fluctuations (Fig. 4.14). Indicator levels should be cautiously interpreted, considering that inverse values of LHE’s were calculated. High values of inverse LHE’s mean low Hölder exponents. As Hölder exponents represent the integrated contribution of long-range dependence and intermittency (equation 4.3), higher 1/LHE’s values (Fig. 4.14, 1993 fishing season) represent lower values of long-range dependence and/or lower values of the intermittency parameter α. For 1996 to 1998 fishing seasons, high values of long-range dependence in TDC residuals quantified by parameter d (Fig. 4.11) contribute to the low observed levels of 1/LHE for those years (Fig. 4.14). Together, this indicator suggests a change in the pattern of long-range dependence and/or intermittency of catches in 1994 and 1995 compared with other years. 127 Figure 4.14. (Upper plot) Time series of total daily catch (TDC) residuals for 1993-1998 fishing seasons. (Lower plot) Inverse Local Hölder Exponent (1/LHE’s, red line) estimated using the KG method is presented for a better visualization of early changes in exponents in relation to biomass fluctuations. Vertical lines separate annual fishing seasons. Blue squares represent vulnerable biomass reconstructed for Pacific Fisheries Management Areas 124 and 125 off the west coast of Vancouver Island (e.g. Fig. 3.1) according to Martell et al. (2000). Shaded area encloses fishing seasons 1994 and 1995, for which major changes in biomass occurred. Biomass levels were localized on the plot at the end of April for each fishing season because fishery independent biomass surveys were conducted at the end of April or beginning of May each year. Note that 1/LHE estimated from TDC tracks well the smooth pink shrimp vulnerable biomass for all years. The 1995 fishing season is in the same direction but displaced about 6 months later, possibly affected by the anomalously high fishing effort (DFO 2001) and catches (Fig. 4. 14) due to high prices for shrimp during 1995 (DFO 2013). The transition of biomass levels from one year to the next is tracked well by 1/LHE for all fishing seasons. This means that it is 128 possible to have an indicator of major changes in the biomass at the beginning and during the progress of each fishing season (Fig. 4.14). If we consider that the dynamics of this stock appear to be driven by recruitment changes, the continuous and sharp decrease in shrimp recruitment (not shown) starting from 1993 up to the 1996 fishing season (Martell et al. 2000) could be the main cause of the decline of the inverse local d indicator for the same time period (Fig. 4.15). From 1996 to 1998, the decline in recruitment continues but at a slower rate (Martell et al. 2000), similar to the trend of the inverse local d indicator (Fig. 4.15). Figure 4.15. (Upper plot) Time series of total daily catch (TDC) residuals for 1993-1998 fishing seasons. (Lower plot) Inverse local d parameter (1/d, red line) with vertical grey lines separating annual fishing seasons. Blue squares represent vulnerable biomass reconstructed for Pacific Fisheries Management Areas 124 and 125 according to Martell et al. (2000). Shaded area encloses fishing seasons 1994 and 1995, for which major changes in biomass occurred. Biomass levels were localized on the plot at the end of April for each fishing season because fishery independent biomass surveys were conducted at the end of April or beginning of May each year. 129 The inverse local d parameter for the early 1993 fishing season fluctuates around 10 and then continuously decreases until the end of the 1995-fishing season when it fluctuates around 3 (Fig. 4.15). This means that TDC exhibited an increasing level of long-range dependence for the 1993-1995 period, stopping its accelerated increase at the end of 1995. For 1996 to 1998 fishing seasons, d reached values slightly greater than 2 (Fig. 4.15), which corresponds to values slightly greater than the upper theoretical limit of the long-range dependence parameter (0 < d < 1/2; Percival and Walden 2000). Using the historical trends of this indicator to detect major changes in vulnerable biomass at the beginning of each fishing season, it is possible to discern between two dominant stages. The first stage exhibits a continuous decrease of vulnerable biomass from 1994 to 1996, and the second stage is characterized by a period of relative stability for 1997 and 1998 fishing seasons (Fig. 4.15) coinciding with a decrease in the exploitation rates (Martell et al. 2000, Walters and Martell 2004), and with the setting of quotas (Rutherford et al. 2004). In contrast, local α values are more variable than LHE’s and local d exponents, and therefore, they were not as useful to discern between different stages of the fishery in terms of biomass levels. In this context, local α only partially tracks the decrease of biomass after the 1995 fishing season, but this indicator does not provide information about the decrease in biomass levels between 1994 and 1995 fishing seasons (Fig. 4.16). According to these results, it appears than the long-range dependence parameter d is controlling the pattern of LHE’s. 130 Figure 4.16. (Upper plot) Time series of total daily catch (TDC) residuals for 1993-1998 fishing seasons. (Lower plot) Local intermittency exponent (α, red line) with vertical grey lines separating annual fishing seasons. Blue squares represent vulnerable biomass reconstructed for Pacific Fisheries Management Areas 124 and 125 according to Martell et al. (2000). Shaded area encloses fishing seasons 1994 and 1995, for which major changes in biomass occurs. Biomass levels were localized on the plot at the end of April for each fishing season because fishery independent biomass surveys were conducted at the end of April or beginning of May each year. 4.3.3 Comparison of early warning indicators The performance of classical statistical indicators proposed to detect early changes in ecosystems was compared with the ability of local Hölder exponents (LHE’s) to detect those changes. Linear patterns were detected for local Hölder exponents (LHE’s), autoregressive lag-1 coefficients (AR1), and spectral ratios (Fig. 4.17a,b). Linear patterns 131 for skewness and kurtosis are not clearly visible, except for an early segment of both time series associated with 1993 and 1994 fishing seasons (Fig. 4.17c,d). Major changes in fishing effort (e.g. Fig. 4.2b) during the 1995 fishing season appear to be well tracked by a change in the growth rate (slope) of local Holder exponents. Figure 4.17. Dynamics of local Hölder exponents (LHE, red line) and fluctuations of total daily catch (TDC) residuals for 1993-1998 fishing seasons (black line) in comparison to a) autoregressive lag-1 coefficient, b) density ratio fluctuations, c) skewness and d) kurtosis (blue lines). Vertical grey lines separate annual fishing seasons. AR(1) coefficients and density ratios at the beginning of this season (Figs. 4.17a,b and 4.18a). Skewness and kurtosis were not sensitive to changes in fishing effort (Fig. 4.17c,d and 4.18b). However, if we focused on the two major biomass declines identified from fishery-independent survey results between 1994-1995 and 1995-1996 reaching biomass levels below the limit reference point (LRP) of 2342.9 ton during the later year (DFO 2012) 132 which pattern is also visible on Figure 4.14, changes in LHE’s estimated from TDC become evident after calculating local slopes (Fig. 4.18a). However, the area under the LHE curve is greater for fishing season 1995 than for fishing season 1994 or 1996, which means that the LHE responds not only to biomass declines but also to the anomalously high fishing effort exerted on the fishery during 1995 (Fig. 4.18a). Figure 4.18. Patterns of change of local slopes for local Hölder exponents (LHE’s, continuous red line) in comparison to a) local slopes of autoregressive lag-1 coefficients (continuous grey line) and local slopes of density ratios (dotted grey line); and b) local slopes of skewness (continuous grey line) and local slopes of kurtosis (dotted grey line) (scaled to fit this figure). Vertical grey lines separate annual fishing seasons. A window equal to 91 observations was used to calculate local slopes. 133 For time periods characterized by stable biomass levels and changing fishing effort (from middle 1997 to end of 1998) (Figs. 4.2b and 4.14), LHE’s showed no major variability or departures from a slope equal to zero (Fig. 4.18), reinforcing the hypothesis that LHE’s can be used to detect changes in biomass in advance of other, more traditional, methods. In contrast, AR(1) coefficients and spectral ratios were highly volatile showing substantial variability during fishing seasons (1997 and 1998) for which total daily catches and biomass were relatively stable (Figs. 4.2c and 4.14). Skewness and kurtosis fluctuate for almost every fishing season (Fig. 4.18b), which invalidates their use to detect early changes in this smooth pink shrimp fishery. 4.4 Discussion 4.4.1 Early detection of fisheries changes The general lack of reference points for Pandalid shrimps poses a risk to the long-term sustainability of shrimp fisheries (Cadrin et al. 2004). Under this framework, empirical indicators can be a valid alternative when traditional fisheries reference points (e.g. spawning stock biomass, SSB) are not available to fisheries managers. The use of empirical indicators has been proven to be useful in detecting early changes of exploited terrestrial and marine ecosystems (Guttal and Jayaprakash 2008, Biggs et al. 2009). As a consequence, they can be considered as valid alternatives for the management of aquatic ecosystems subjected to different fish predation levels and harvest rates (Biggs et al. 2009, Litzow et al. 2013). The importance of the early detection of major changes underlying a fishery has increased over time. Considering that factors driving the recovery of overexploited fish stocks are 134 multiple and uncertain, and that most of these stocks remain below target biomass levels (Neubauer et al. 2013), early detection of major changes in these stocks that are occurring now may be more valuable than recovery of these fish stocks at some later date from an economic and ecological standpoint. However, prediction of marine fisheries dynamics is extremely limited due to the complex nonlinear dynamics that emerge from the coupling of social and ecological systems (Glaser et al. 2014). In this chapter the focus has been on efforts to develop self-similar indicators in order to provide early detection of major changes underlying the smooth pink shrimp fishery off the west coast of Vancouver Island, Canada, as a management tool to help with this, and potentially other, fisheries. One main advantage of our proposed 1/LHE indicator is that it can be used for “real time” prediction. For example, to predict the 1994/1995 vulnerable biomass change the historical catch data before 1994/1995 should be available, but not data after the change. In addition, regular updates to predictions can be conducted with the latest available catch data. We are aware that generic signals warning of changes in ecosystems, including exploited marine ecosystems, may not exist (Boettiger and Hastings 2012). Different exploited ecosystems, however, may still exhibit some common properties that can be used to construct indicators to detect impending major changes. In this context, the smooth pink shrimp fishery from the west coast off Vancouver Island and the southern hake fishery off Chile (Montes 2004) exhibit changes in self-similar patterns of variability (long-range dependence and intermittency) before major changes in the fishery take place. Intermittency, measured as an increasing pattern of variability, and skewness have been used to detect in advance the collapse of the crustacean fishery of the Gulf of Alaska and 135 Bering Sea (Litzow et al. 2013). However, use of these indicators by fisheries managers should wait until the empirical validation of their utility on each specific ecosystem. It is important to note that after comparing LHE’s obtained for smooth pink shrimp TDC time series (which varies between approximately 0.82 and 0.92, Fig. 4.10) with those values obtained for the simulated time series (Fig 4.5 to Fig. 4.9), the estimation of H using the KG method can be considered unbiased. This means that in simulated cases for which H is between the ranges of 0.82-0.92 (e.g. d=0.2679 in Fig. 4.5 or d=0.2941 in Fig. 4.6), the mean value of H represented on the box-plot using the KG method was very close to the theoretical reference line. For the full range of LHE’s obtained from the smooth pink shrimp fishery we can be confident that this indicator is unbiased independent of the trends presented by the fishery, its degree of variability or the percentage of missing observations; estimated LHE values are therefore comparable through the history of the fishery and represent major changes in this fishery. For the pink shrimp fishery off the west coast of Vancouver Island, changes in biomass levels and fishing effort are tracked by changes in LHE’s according to equation (4.3). We believe that, as major changes in biomass levels occur at a slower rate involving multiple time scales than changes in fishing effort, biomass changes are better tracked by changes in the LRD parameter d than by the intermittency parameter α. On the contrary, changes in fishing effort were detected better by the intermittency parameter α, which quantifies fast departures from the mean level of effort. This conclusion is supported by our results, because when only the local d parameter was used to detect early changes, major changes in biomass levels are well tracked by the LRD parameter d (Fig. 4.15). When only the intermittency parameter α was used to track changes of biomass, this indicator was 136 substantially affected by the high fishing effort exerted on the fishery during 1994 and 1995 fishing seasons achieving its theoretical upper limit (α=2) for these years (Fig. 4.16). 4.4.2 Long-range correlations in fish recruitment and local Hölder exponents (LHE’s) Recruitment time series of different species of shrimps and from other invertebrates show self-similar patterns of variability, which means that year-to-year recruitment shows long-range correlations (Niwa 2006). Good and bad recruitment years are clustered in time, and their effects do not disappear from one year to the next. This clustered effect in recruitment could be the reason why smooth pink shrimp total daily catches are correlated at inter-annual time scales, pushing the Hurst coefficient estimated from catches towards values greater than 0.5 (i.e. different from a white noise, or uncorrelated, time series; Beran et al. 2013). It appears that the stock dynamics of smooth pink shrimps off the west coast of Vancouver Island is driven by recruitment (Martell et al. 2000). Considering that highly variable yields of pink shrimp populations off the west coast of Vancouver Island and off the Oregon/Washington coast are strongly dependent on annual recruitment (1-year old shrimps; Martell et al. 2000, Hannah 2011), time series of catches could be considered, to some extent, representative of recruitment success. Smooth pink shrimp recruitment time series off the west coast of Vancouver Island (WCVI) showed a period of relative stability at a high level over the 1987 to 1994 fishing seasons, then suffered a substantial decline between 1994 and 1995, and continued to decrease until 1997 at a slower rate (Martell et al. 2000). Inverse local Hölder exponents (1/LHE’s) calculated using total daily catches show dynamics similar to the recruitment 137 anomaly time series (Martell et al. 2000) and also to the vulnerable biomass (Fig. 4.13) which supports the hypothesis that 1/LHE’s can be used to estimate recruitment success before the start of the fishing season using commercial catches. After showing sustained juvenile survival rates (positive recruitment anomaly) from 1991, negative recruitment anomalies were estimated for 1994 in Pacific Fisheries Management Areas (PFMA) 124 and 125 (Martell 2002), which coincided with the sharp decline in 1/LHE’s between 1994 and 1995 (Fig. 4.19) and was also detected by the local slopes plot (Fig. 4.11). Early changes of recruitment regimes like this using 1/LHE’s should be confirmed with the analysis of longer time series showing multiple alternative states with positive and negative recruitment anomalies. 138 Figure 4.19 (Upper plot) Inverse local Hölder exponents (1/LHE, red line) and vulnerable biomass for PFMA 124 and 125 (Martell et al. 2000). (Middle plot) Recruitment anomalies estimated for PFMA 124. (Lower plot) Recruitment anomalies estimated for PFMA 125 (Martell 2002). 4.4.3 Comparison of early warning indicators Many indicators have been constructed based on the critical slowing down (CSD) concept (Strogatz 2001) which can be quantified statistically as a rise in variance and autocorrelation, or as a shift in skewness or kurtosis before a major ecosystem change takes place (van Nes and Scheffer 2007, Brock et al. 2008, Guttal and Jayaprakash 2008). Quantitatively, no indicators had outperformed others (Dakos et al. 2012, Boettiger et al. 2013). We believe that our self-similar LHE indicator has the capacity to estimate both aspects of time series dynamics, long-range autocorrelation and variability, thereby 139 integrating both into a single index (Taqqu and Teverovsky 1998, Frankze et al. 2012). For that reason, the LHE indicator can better track changes in exploited ecosystems than previously proposed indicators. We want to remark that self-similar indicators can be considered as a broad class of indicators, which includes statistical indicators, such as kurtosis and skewness, particular cases for specific values of α intermittency levels. At the same time, AR(1) estimates the correlation between observations at lag=1 among multiple possible estimates of correlation if larger time lags are included. Autocorrelation coefficient at lag-1, the most commonly used time series property, is a measure of short-term correlation between observations (Beran et al. 2013), but ignores dependence between distant observations, i.e. long-range dependence (LRD). This latter property offers a greater possibility to provide early detection of ecosystem changes than using lag-1 autocorrelation derived indicators. This is because as ecosystem properties start to change, the effect they have on the first observations subjected to this change propagates through time modifying the degree of LRD or the extent of scaling ranges. Obviously, the greater the scaling range detected in normal or baseline conditions, the greater the possibility they offer to detect in advance ecosystem changes. Most indicators fail in detecting ecosystem changes in advance if the change in drivers is fast in comparison with the dynamics of the ecosystem components. For example, if harvest rates applied in aquatic ecosystems increase too fast, indicators fail to detect such changes before they actually occur (Carpenter et al. 2008, Biggs et al. 2009). In our case, for the smooth pink shrimp TDC time series, abrupt changes in fishing effort during the 1995 fishing season (Fig. 4.2b), which appear to be not totally filtered out, partially altered the 140 dynamics of 1/LHE and 1/d indicators. Both indicators show a period of relative stability before continuing their general trend of decline (Figs. 4.14 and 4.15, respectively). The analyses of smooth pink shrimp catch rates within each fishing season may give information about changes in stock size as the fishery proceeds (Martell et al. 2000). For other invertebrate fisheries (Roa-Ureta and Arkhipkin 2007, Murray et al. 2013) daily catch rates at the beginning of the fishing season better represent biomass than catches rates at the middle or at the end of the fishing season. In this context, self-similar indicators can be used to help in the estimation of reliable initial quotas for the smooth pink shrimp fishery for any fishing season (starting in April of year t) based on the values of 1/LHE obtained for TDC available up to October of year t-1, which means five months in advance of the start of the season. It is important to remember that daily smooth pink shrimp TDC time series per se are not an index of recruitment and/or biomass. However, the dependence between historical observations of total daily catches, extracted at particular time scales over a long-term horizon, can be used to construct an index capable of providing early detection of changes underlying this fishery. The use of the commercial fisheries data enable construction of indicators that permit the detection of changes on a time scale that is useful in, and can be used for, the implementation of in-season smooth pink shrimp fisheries management strategies. Considering that the type of information used in the elaboration of self-similar-based indicators is available for many other fisheries, and most indicators usually fail to detect ecosystem changes when annual data are used to predict such changes (Perreti and Munch 2012), fine scale catch data for other fisheries should be used to test the conditions, performance and generality of proposed early-warning indicators. 141 There has been much effort in recent years on detecting large-scale abrupt changes in ecosystems and critical thresholds over time. The ultimate goal of this line of research is give managers the chance to intervene in order to avoid the collapse of ecosystems and fish populations. Several early warning indicators have been used in the detection of major changes in aquatic ecosystems using single statistical time series properties that change when ecosystems get closer to major perturbations. Based on results obtained in Chapter 2 and Chapter 3, we opened to the possibility that multifractality can vary through time, which was confirmed in Chapter 4. Changes in multifractal dynamics of total daily catches across fishing seasons is a property of this fishery that was detected and quantified in Chapter 4, when a time series indicator able to detect long-range dependence (LRD), intermittent variability and their changes through time, was constructed. One of our proposed indicators, the inverse of the Local Hölder Exponent (1/LHE), integrates LRD and intermittency into a single coefficient, allowing the detection of changes in both statistical time series properties simultaneously in a localized manner. This indicator has proven to be robust to the presence of missing observations and non-Gaussian time series fluctuations, and has the capacity to track well the smooth pink shrimp vulnerable biomass for the study period. This result means that it is possible to have an indicator of major changes in biomass at the beginning and during the progress of each fishing season using only time series of catches. Compared to four statistical indicators (autoregressive lag-1 coefficient, spectral ratios, skewness, kurtosis) used to detect early changes in ecosystems, our indicator shows the capacity to track changes better in smooth pink shrimp biomass than previously proposed indicators, responding also to the anomalously high fishing effort exerted on the fishery. 142 Chapter 5: Summary and Conclusions In this thesis, the temporal dynamics of smooth pink shrimp catches off the west coast of Vancouver Island and oceanographic time series from the same area were analyzed from several perspectives under the framework of fractal theory. The main objectives of this thesis were: 1) To determine and quantify fractal patterns of variability of smooth pink shrimp total daily catches off the west coast of Vancouver Island. 2) To determine and quantify fractal patterns of variability of oceanographic time series. 3) To determine and quantify bivariate fractal or multifractal relationships between oceanographic and smooth pink shrimp fisheries time series. 4) To evaluate the usefulness of fractal parameters as indicators of change in the dynamics of Pandalus jordani fishery. To achieve these objectives, a sound understanding of fractal theory was necessary. In consequence, we started on Chapter 1 with an introduction to fractals, including a detailed description of its geometrical properties, known as self-similarity and self-affinity. Then, the concept of fractal dimension (D), as a measure of structural complexity, was introduced. The fundamental power-law relationship between the length of a fractal object and the scale unit used to measure it was explained using equation (1.1). This latter equation represents the theoretical basis of fractal geometry, and explains the way that fractal dimension can be estimated. Then, the self-similar property underlying equation (1.1) was used to extend the use of fractals from the spatial to the temporal domain, which is the focus of this thesis. 143 Fractional Gaussian noise (fGn) and fractional Brownian motions (fBm) were introduced in order to understand their statistical properties. In addition, these latter functions were used to explain how the Hurst coefficient (H) regulates the correlation structure between distant observations, and how it quantifies self-affinity in long-range dependent time series. This coefficient was estimated, using different methods, in fisheries and oceanographic time series in Chapter 2, Chapter 3 and Chapter 4. Special attention was given to the self-affine property of these latter functions, which means that their fractal properties are expressed on average conditions. For this reason, we focused on the analysis of stochastic fractal time series using an average quantity estimated over several time lags to obtain an estimation of the Hurst coefficient (H). When necessary, fishing day’s time series were created after the stitching of smooth pink shrimp catches of adjacent fishing seasons to create continuous time series. The time lag used in this thesis was one fishing day and represents the scale unit and time lag expressed in equation (1.8). Fractional Gaussian noise (fGn) and fractional Brownian motion (fBm) functions with increasing Hurst coefficients were simulated to illustrate the difference between no long-range dependence (white noise, H=0.5) and pronounced long-range dependence (H=0.9) for stationary (fGn) and non-stationary (fBm) time series. In the second half of Chapter 1, multifractals were introduced and explained as a special type of self-affine time series with a hierarchy of subsets, each having a different fractal dimension (Seuront 2010). Then, the order of singularity concept was introduced, representing a defined intermittent level quantified by parameter C1 using the Universal Multifractal Model (UMM, Lovejoy et al. 2001; equation 1.9). This latter parameter, in conjunction with the multifractal parameter α, also known as the Levy index, regulates the 144 frequency of jumps or intermittency level of a multifractal time series. Both parameters were useful, in conjunction with the Hurst coefficient, for the detection of early changes in the temporal dynamics of smooth pink shrimp catches conducted in Chapter 4. At the end of Chapter 1, examples of applications of fractals in benthic and planktonic populations were presented. This is not an exhaustive review, but has the intention to highlight the great advances conducted in fractal research in marine spatial ecology in comparison to studies focused on the temporal domain. At the end of this chapter, five case studies were presented which focused on quantifying the long-range dependence (LRD) temporal dynamics of diverse fish stocks. Objective 1 was accomplished in Chapter 2. Multifractal patterns of variability were detected and quantified in the chapter using three continuous smooth pink shrimp fishing seasons (1994, 1995 and 1996) from the west coast off Vancouver Island, Canada. Long-range dependence (LRD) and intermittent dynamics were quantified using the Universal Multifractal Model (UMM, Lovejoy et al. 2001) and the Structure Function methodology (Seuront et al. 1999). A scaling range between 16 and 120 fishing days was selected using an automatic adaptive scaling range selection algorithm (Xia et al. 2005), for which, a Hurst coefficient H=0.8778 was estimated. A comparison between a fractal and a multifractal model was conducted using a likelihood ratio test for nested models (Burnham and Anderson 2002). The multifractal model was selected as significantly better than the monofractal model. This result means that long-range dependence was detected for smooth pink shrimp catches, but catches are more complex than can be described by only a single long-range dependent parameter (H). In other words, intermittency behaviour quantified by 145 parameter C1=0.0429 and α=1.5073, explained the occurrence of rare and infrequent observations and the non-Gaussian probability distribution of catches, respectively. The usefulness of the Universal Multifractal Model was highlighted at the end of this chapter from a fisheries management perspective. This model can be used to predict the probability of achieving different thresholds of catches, and to interpret the occurrence of occasional large catches. Furthermore, it is possible to use this model to forecast smooth pink shrimp catch values at intra-annual time scales, which helps in predicting when Total Allowable Catches will be reached and fishing areas may be closed. In-season estimation of future variability of catches should be considered as an additional tool to help in managing this fishery, especially after the completion of fishery-independent biomass surveys and the setting of preliminary catch ceilings. Objective 2 was accomplished on Chapter 3 when fractal patterns were quantified in sea surface temperature (SST) and wind stress (WS) time series from the west coast off Vancouver Island. Using Detrended Fluctuation Analysis (DFA, Kantelhardt et al. 2006), long-range dependence was estimated for sea surface temperature (H=0.9675) for an extended time series that covers from 1940 to 2013; and for wind stress (H=0.7538) for the 1958-1998 time series. In this chapter, the variability of Hurst coefficients across 70 years was quantified, and its relation to anomalous oceanographic conditions was discussed. These results are in agreement with those found for the wider North-East Pacific Ocean, which reinforces the potential predictability skill that SST has over long-term horizons outperforming short-range auto-regressive memory models (Zhu et al. 2010). 146 Working under the hypothesis that long-range dependence will be preserved through the history of the fishery, a longer time series of catches which covers 12 stitched fishing seasons (1987-1998), in comparison to time series used in Chapter 2 (3 stitched fishing seasons, 1994-1996), was analyzed. In this context, long-range dependence was quantified using Detrended Fluctuation Analysis (DFA), which allowed the estimation of long-range dependence in total daily catch time series for a time period contaminated with data gaps and different trends. Two scaling ranges were obtained. The first scaling range extends between 16 and 120 fishing days, and the second between 120 and 642 days. The first scaling range coincides with that obtained when H was estimated for the three stitched fishing season time series of catches in Chapter 2 (H=0.8779, N=642). The Hurst coefficient obtained in Chapter 3 for this longer series is slightly lower (H=0.8568, N=2568), corroborating that long-range dependence is present in smooth pink shrimp total daily catches. However, weak evidence of multifractality was detected in Chapter 3 in the contrary to results obtained in Chapter 2. Thus in Chapter 4, changes in multifractal dynamics across the history of the fishery were assessed. In the second half of Chapter 3, Objective 3 was accomplished. Long-range cross-correlations between total daily catch anomaly and filtered sea surface temperature time series (ho,f=0.8580) and between total daily catch anomaly and filtered wind stress time series (ho,f=0.8315) were detected and quantified. Much effort was conducted in this chapter on removing periodic oscillations in oceanographic and fisheries time series, because these deterministic trends seriously distort the estimation of LRD parameters (Montanari et al. 1999). The continuous wavelet transform (CWT, Torrence and Compo 1998) was used to remove oceanographic time series variability associated with the 147 distorting effect of quasi-biennial oscillations (2-3 years) and El Niño Southern Oscillation (5-7 years) that have been well described for this ecosystem (Ware 1995). This latter step is just as important as the detection of long-range dependence, and it must be conducted carefully before fractal parameters are estimated. In Chapter 4, the last objective of this thesis (Objective 4) was met when localized indicators defined as local Hölder exponents (LHE’s) were constructed, and their usefulness in tracking in-season dynamics of smooth pink shrimp total daily catches (TDC) was tested. Time series with different degrees of long-range dependence (d), intermittency levels (α), and missing observations were simulated using a FARIMA model (Stoev and Taqqu 2004) to test the uncertainty and bias of two methods used for estimating H. The KG method examined in this research for this purpose performs better than the QQ-ML estimator. This result gave us confidence when LHE’s were calculated using the KG method because these exponents are unbiased, independent of the trends presented by the smooth pink shrimp fishery and not affected by the variability of catches or the percentage of missing observations. Local Hölder exponents (LHE’s) integrate the contribution of long-range dependence and intermittency into a single index, providing the ability to detect early changes underlying this fishery. This was demonstrated using historical biomass as an indicator of the status of shrimp populations. It appears that multifractality varies across the history of the smooth pink fishery, therefore multifractal dynamics should be quantified in a localized manner. However, inter-annual changes in shrimp biomass from one fishing season to the next were well tracked by continuous changes in the fractal-based indicator (1/LHE) constructed in Chapter 4. The usefulness of indicators specified in Objective 4 was demonstrated in this 148 Chapter when the index proposed here (1/LHE) estimated from total daily catches shows its capacity to track well the smooth pink shrimp vulnerable biomass for all years (1993-1998). In consequence, 1/LHE index was proposed to be used as an indicator of change in shrimp biomass at the beginning and during the progress of each fishing season. The results of this chapter were compared with statistical indicators used to detect in advance changes in exploited ecosystems (Dakos et al. 2012). Two useful properties were detected for the 1/LHE index: it reacts faster to changes in biomass than the other statistical indicators examined here, and it is relatively insensitive during periods characterized by stable biomass levels and changing fishing effort. In summary, 1/LHE is proposed as a valid management tool to help with this, or potentially other, fisheries. According to the results obtained in Chapter 4, pink shrimp catches appears to be more consistent with the dynamics of multifractional processes in comparison to a multifractal process. For a multifractional process, a family of different exponents are used to estimate how self-similarity varies through time (Sheng et al. 2012a,b). Further studies should be conducted to prove if 1/LHE index tracks well changes in biomass of other exploited fish populations from the Northeast Pacific and other marine ecosystems before claiming its generality. The ability of this indicator to detect marked inter-annual changes underlying the fishery has been proven; however, the early detection of major historical collapses suffered by well-documented fisheries such as the Atlantic Cod fishery (Myers et al. 1997) should be one of the next steps in testing the general applicability of this indicator. In addition, the response of 1/LHE indicator to changes in harvest technology and market forces should be quantified, as it was already conducted for other statistical early warning indicators (Richter and Dakos 2015). 149 Usually, catch data is regularly collected by fisheries agencies, and taking into account the valuable information provided by 1/LHE indicator using daily catches, the collection of catch data on this fine scale is strongly encouraged for this and other fisheries. We are aware that today the use of daily logbook data for the elaboration of this indicator could be a major limitation for many fisheries. In consequence, valuable information could be obtained if future work is conducted focusing on quantifying the effect of using coarser temporal scales (weekly, monthly data) in the construction of the 1/LHE indicator. The construction of 1/LHE is straightforward, however, it should not be applied in a blinded manner. The detection and removal of deterministic trends (periodic oscillations) should be removed carefully before the application of this indicator. The early detection of threshold levels over which a fishery moves to an undesirable or deplete state should be analyzed under the framework of the proposed 1/LHE indicator. In this context, the quantification of time available to fishery managers to intervene, according to the signal given by the 1/LHE indicator, should be quantified in future studies to prove the utility of this indicator. 150 References Abry, P., Gonçalves, P. and Véhel, J.L. 2009. Scaling, Fractals and Wavelets. Wiley-ISTE Ltd, London, UK. Addison, N. 2002. The Illustrated Wavelet Transform Handbook. Introductory Theory and Applications in Science, Engineering, Medicine and Finance. Institute of Physics Publishing, London, 400 pp. Ashkenazy, Y. and Gildor, H. 2009. Long-range temporal correlations of ocean surface currents. Journal of Geophysical Research 114: C09009. Akçakaya, H.R., Halley, J.M. and Inchausti, P. 2003. Population-level mechanisms for reddened spectra in ecological time series. Journal of Animal Ecology 72: 698-702. Amos, C.L., Martino, S., Sutherland, T.F. and Rashidi, T.A. 2015. Sea surface temperature trends in the coastal zone of British Columbia, Canada. Journal of Coastal Research 31: 434-446. Avnir, D., Biham, O., Lidar, D. and Malcai, O. 1998. Is the Geometry of Nature Fractal? Science 279: 39-40. Azovsky, A.I., Chertoprood, M.V., Kucheruk, N.V., Rybnikov, P.V. and F.V. Sapozhnikov. 2000. Fractal properties of spatial distribution of intertidal benthic communities. Marine Biology (136): 581-590. Baumgartner, T.R., Soutar, A. and V. Ferreira-Bartrina. 1992. Reconstruction of the history of Pacific sardine and northern anchovy populations over the past two millennia 151 from sediments of the Santa Barbara basin, California. California Cooperative Oceanic Fishery Investment Report (33): 24-40. Barbosa, S.M. 2011. Testing for deterministic trends in global sea surface temperature. Journal of Climate 24: 2516-2522. Barbosa, S.M., Silva, M.E. and Fernandes, M.J. 2008. Time series analysis of sea level records: characterising long-term variability. In Donner R.V. and S.M. Barbosa (eds.) Nonlinear Time Series Analysis in the Geosciences. Applications in Climatology, Geodynamics and Solar-Terrestrial Physics, Springer-Verlag, Berlin, pp. 157-173. Bailey, K.M., Hollowed, A.B. and Wooster, W.S. 2004. Complexity of marine fisheries dynamics and climate interactions in the Northeast Pacific Ocean. In Stenseth N. C., Ottersen, G., Hurrel, J.W. and A. Belgrano (Eds.), Marine Ecosystems and Climate Variation, Oxford University Press, Oxford, UK, pp. 147-151. Batt, R.D., Carpenter, S.R., Coleb, J.J., Pace, M.L. and Johnson, R.A. 2013. Changes in ecosystem resilience detected in automated measures of ecosystem metabolisms during a whole-lake manipulation. Proceedings of the National Academy of Sciences 110: 17398-17403. Beamish, R.J., McFarlane, G.A. and King, J.R. 2005. Migratory patterns of pelagic fishes and possible linkages between open ocean and coastal ecosystems off the Pacific coast of North America. Deep-Sea Research II 52: 739-755. 152 Beran, J., Feng, Y., Ghosh, S. and Kulik, R. 2013. Long-Memory Processes. Probabilistic Properties and Statistical Methods. Springer-Verlag, Berlin, 884 pp. Bertrand, S., Burgos, J.M., Gerlotto, F. and J. Atiquipa. 2005. Lévy trajectories of Peruvian purse-seiners as an indicator of the spatial distribution of anchovy (Engraulis ringens). ICES Journal of Marine Science (62): 477-482. Biggs, R., Carpenter, S.R. and Brock, W.A. 2009. Turning back from the brink: Detecting and impending regime shift in time to avert it. Proceedings of the National Academy of Sciences 106: 826-831. Biham, O., Malcai, O., Lidar, D. and Avnir, D. 1998a. Is nature fractal. Science 279: 785-786. Biham, O., Malcai, O., Lidar, D. and Avnir, D. 1998b. Fractality in Nature. Science 279: 1615-1616. Bjørnstad, O., Nisbet, R.M. and Fromentin, J.M. 2004. Trends and cohort resonant effects in age-structured populations. Journal of Animal Ecology 73: 1157-1167. Boettiger, C. and Hastings, A. 2012. Early warning signals and the prosecutor’s fallacy. Proceedings of the Royal Society of London B 279: 4734-4739. Boettiger, C., Ross, N. and Hastings, A. 2013. Early warning signals: the charted and uncharted territories. Theoretical Ecology 6: 255-264. Boucharel, J., Dewitte, B., Garel, B. and Penhoat, Y. 2009. ENSO’s non-stationary and non-Gaussian character: the role of climate shifts. Nonlinear Processes in Geophysics 16: 453-473. 153 Boufadel, M.C., Lu, S., Molz, F.J. and Lavallée, D. 2000. Multifractal scaling of the intrinsic permeability. Water Resources Research 36: 3211-3222. Boutillier, J.A., Perry, R.I., Waddell, B., and Bond, J. 1997. Assessment of the offshore Pandalus jordani trawl fishery off the west coast of Vancouver Island. Pacific Stock Assessment Review Committee (PSARC) Working Paper No. I-97-11. Brock, W.A., Carpenter, S.R. and Scheffer, M. 2008. Regime shifts, environmental signals, uncertainty and policy choice. In: Norberg, J. and G. Cumming (eds). A theoretical framework for analyzing socio-ecological systems, Columbia, New York, USA, pp. 180-206. Buckley, T.W. and Livingston, P.A. 1997. Geographic variation in the diet of Pacific hake, with a note on cannibalism. California Cooperative Oceanic Fisheries Investigations Reports 38: 53-62. Burnham, K.P. and Anderson, D.R. 2002. Model Selection and Multimodel Inference. A Practical Information-Theoretic Approach. Second Edition, Springer-Verlag, New York. Butler, T.H. 1980. Shrimps of the Pacific Coast of Canada. Canadian Bulletin of Fisheries and Aquatic Sciences 202. Bradbury, R.H., Reichelt, R.E. and D.G. Green. 1984. Fractals in ecology: methods and interpretation. Marine Ecology Progress Series (14): 295-296. 154 Bunde, A., Havlin, S., Kantelhardt, J.W., Penzel, T., Peter, J.-H. and K. Voigt. 2000. Correlated and Uncorrelated Regions in Hearth-Rate Fluctuations during Sleep. Physical Review Letters (85): 3736-3739. Burrough, P.A. 1981. Fractal dimensions of landscape and other environmental data. Nature 294: 240-242. Bylhouwer, B., Ianson, D. and Kohfeld, K. 2013. Changes in on the onset and intensity of wind-driven upwelling and downwelling along the North American Pacific coast. Journal of Geophysical Research 118: 2565-2580. Caccia, D.C., Percival, D., Cannon, M.J., Raymond, G. and Bassingthwaighte, J.B. 1997. Analyzing exact fractal time series: evaluating dispersional and rescaled range. Physica A 246: 609-632. Cadrin, S.X., Boutillier, J.A. and Idoine, J.S. 2004. A hierarchical approach to determining reference points for Pandalid shrimp. Canadian Journal of Fisheries and Aquatic Sciences 61: 1373-1391. Cannon, M.J., Percival, D.B., Caccia, D.C., Raymond, G.M. and Bassingthwaighte, J.B. 1997. Evaluating scaled windowed variance methods for estimating the Hurst coefficient of time series. Physica A 241: 606-626. Carpenter, S.R. and Brock, W.A. 2006. Rising variance: a leading indicator of ecological transition. Ecology Letters 9: 311-318. Carpenter, S.R., Brock, W.A., Cole, J.J., Kitchell, J.F. and Pace, M.L. 2008. Leading indicators of trophic cascades. Ecology Letters 11: 128-138. 155 Carpenter, S.R., Cole, J.J., Pace, M.L., Batt, R., Brock, W.A., Cline, T., Coloso, J., Hodgson, J.R., Kitchell, J.F., Seekell, D.A., Smith, L. and Weidel, B. 2011. Early Warnings of Regime Shifts: A Whole-Ecosystem Experiment. Science 332: 1079-1082. Carpenter, S.R., Brock, W.A., Cole, J.J. and M.L. Pace. 2014. A new approach for rapid detection of nearby thresholds in ecosystem time series. Oikos 123: 290-297. Cazelles, B., Chavez, M., Berteaux, D., Ménard, F., Vik, J.O., Jenouvrier, S and Stenseth, N. Chr. 2008. Wavelet analysis of ecological time series. Oecologia 156: 287-304. Chatfield, Ch. 2000. The Analysis of Time Series. An Introduction. Chapman & Hall/CRC, 283 pp. Chamoli, A., Bansal, A.R. and Dimri, V.P. 2007. Wavelet and rescaled range approach for the Hurst coefficient for short and long time series. Computers & Geosciences 33: 83-93. Chen, Z., Ivanov, P.C., Hu, K. and Stanley, H.E. 2002. Effect of nonstationarities on detrended fluctuation analysis. Physical Review E 65: 041107. Chen, Z., Hu, K., Carpena, P., Bernaola-Galvan, P., Stanley, H.E. and Ivanov, P. Ch. 2005. Effect of nonlinear filters on detrended fluctuation analysis. Physical Review E 71: 011104: 1-11. Crawley, M.J. 2007. The R Book. John Wiley & Sons, Chichester. Cryer, J.D. and K.-S. Chan. 2008. Time Series Analysis with Applications in R. Second Edition, Springer Science+Business Media, LCC, New York, NY, USA. 156 Cummins, P.F. and Masson, D. 2014. Climatic variability and trends in the surface waters of coastal British Columbia. Progress in Oceanography 120: 279-290. Dagendorf, S., Rybski, D., Mudersbach, C., Müller, A., Kaufmann, E., Zorita, E. and Jensen, J. 2014. Evidence for long-term memory in sea level. Geophysical Research Letters 41: 5530-5537. Dakos, V., Scheffer, M., van Nes, E.H., Brovkin, V., Petoukhov, V. and Held, H. 2008. Slowing down as an early signal for abrupt climate change. Proceedings of the National Academy of Sciences 105: 14308-14312. Dakos, V., Kéfi, S., Rietkerk, M., van Nes, E.H. and Scheffer, M. 2011. Slowing down in spatially patterned ecosystems at the brink of collapse. The American Naturalist 177: E153-E166. Dakos, V., Carpenter, S.R., Brock, W.A, Ellison, A.M., Guttal, V., Ives, A.R., Kéfi, S., Livina, V., Seekell, D.A., van Nes, E.H. and Scheffer, M. 2012. Methods for detecting early warnings of critical transitions in time series illustrated using simulated ecological data. Plos One 7: e41010. Dakos, V., van Nes, E.H. and Scheffer, M. 2013. Flickering as an early warning signal. Theoretical Ecology 6: 309-317. Davis, A., Marshak, A., Wiscombe, W. and Cahalan, R. 1994. Multifractal characterization of nonstationarity and intermittency in geophysical fields: Observed, retrieved, or simulated. Journal of Geophysical Research 99 (D4): 8055-8072. 157 Davis, A., Marshak, A., Wiscombe, W. and Cahalan, R. 1996. Multifractal characterizations of intermittency in nonstationary geophysical signals and fields. A model based perspective on ergodicity issues illustrated with cloud data. In Currents Topics in Nonstationary Analysis. Edited by G. Treviño, J. Hardin, B. Douglas and E. Andreas. World-Scientific, Singapore, pp. 97-158. Davies, R.B. and Harte, D.S. 1987. Test for Hurst effect. Biometrika 74: 95-101. Delignieres, D., Ramdani, S., Lemoine, L., Torre, K., Fortes, M. and Ninot, G. 2006. Fractal analyses for ‘short’ time series: A re-assessment of classical methods. Journal of Mathematical Psychology 50: 525-544. Denny, M.W., Helmuth, B., Leonard, G.H., Harley, Ch.D.G., Hunt, L.J.H. and Nelson, E.K. 2004. Quantifying scale in ecology: lessons from a wave-swept shore. Ecological Monographs 74: 513-532. Deriche, J.A. and Tewfik, A.H. 1993. Maximum likelihood estimation of the parameters of discrete fractionally differenced Gaussian noise process. IEEE Transactions on Signal Processing 41: 2977-2989. DeVries, T.J. and Pearcy, W.G. 1982. Fish debris in sediments of the upwelling zone off central Peru: a late Quaternary record. Deep-Sea Research 28: 87-109. DFO, 2001. Fish Stocks of the Pacific Coast. Fisheries and Oceans Canada, Canada. DFO, 2011. Pacific Region, Integrated Fisheries Management Plan, Shrimp by Trawl, April 1, 2011 to March 31, 2012. http://www.dfo-mpo.gc.ca/Library/343050.pdf. (accessed 26 September 2011). 158 DFO, 2012. Shrimp Survey Bulletin 12-01, 11 pp. http://www.pac.dfo-mpo.gc.ca/fm-gp/commercial/shellfish-mollusques/shrimp-pcrevette/index-eng.html. (accessed 20 November 2014). DFO, 2013. Pacific Region, Shrimp Trawl Integrated Fisheries Management Plan 2013-2014. Fisheries and Oceans Canada, Canada, 129 pp. (accessed April 15 2015). Dijkstra, H.A. 2013. Analysing data from stochastic dynamical systems. In: Nonlinear Climate Dynamics, Cambridge University Press, pp. 83-115. Dolan, S.S., Bean, C.J. and Riollet, B. 1998. The broad-band fractal nature of heterogeneity in the upper crust from petrophysical logs. Geophysical Journal International 132: 489-507. Dorn, M.W. 2001. Fishing behavior of factory trawlers: a hierarchical model of information processing and decision-making. ICES Journal of Marine Science 58: 238-252. Downton, M.W., and Miller, K.A. 1998. Relationships between Alaskan salmon catch and North Pacific climate on interannual and interdecadal time scales. Canadian Journal of Fisheries and Aquatic Sciences 55: 2255-2265. Eales, J. and Wilen, J.E. 1986. An examination of fishing location choice in the pink shrimp fishery. Marine Resource Economics 2: 331-351. Efron, B. and Tibshirani, R.J. 1993. An Introduction to the Bootstrap. Monographs on Statistics and Applied Probability 57, Chapman & Hall/CRC, Boca Raton, Florida. 159 Eke, A., Herman, P., Kocsis, L. and Kozak, L.R. 2002. Fractal characterization of complexity in temporal physiological signals. Physiological Measurement 23: R1-R38. Evertsz, C.J.G. and Mandelbrot, B.B. 1992. Multifractal Measures. In: Peitgen, H.-O., Jürgens, H. and D. Saupe (Eds). Chaos and Fractals: New Frontiers of Science. Springer-Verlag, New York, pp. 921-953. Fang, W. and Hsieh, W.W. 1993. Summer sea surface temperature variability off Vancouver Island from satellite data. Journal of Geophysical Research 98 (C8): 14391-14400. Fattahi, M.H., Talebbeydokhti, N., Rakhshandehroo, G.R., Shamsai, A. and Nikooee, E. 2011. The robust fractal analysis of time series: concerning signal class and data length. Fractals 19: 29-49. Faucher, M., Burrows, W.R. and Pandolfo, L. 1999. Empirical-statistical reconstruction of surface marine winds along the western coast of Canada. Climate Research 11: 173-190. Faÿ, G., Moulines, E., Roueff, F. and Taqqu, M. 2009. Estimators of long-memory: Fourier versus wavelets. Journal of Econometrics 151: 159-177. Feder, J. 1988. Fractals. Plenum Press, New York. Ferriere, R. and Cazelles, B. 1999. Universal power laws govern intermittent rarity in communities of interacting species. Ecology 80: 1505-1521. 160 Ficker, T., and Benešovský, P. 2002. Deterministic fractals. European Journal of Physics 23: 403-408. Finn, D., Lamb, B., Leclerc, M.Y., Lovejoy, S., Pecknold, S. and Schertzer, D. 2001. Multifractal analysis of line-source plume concentration fluctuations in surface-layer flows. Journal of Applied Meteorology 40: 229-245. Fisher, K.E., Wiebe, P.H. and Malamud, B.D. 2004. Fractal characterization of local hydrographic and biological scales of patchiness on Georges Bank. In Handbook of Scaling Methods in Aquatic Ecology. Measurements, Analysis, Simulation. Edited by L. Seuront and P.G. Strutton. CRC Press, Boca Raton, Florida, pp. 297-319. Foreman, M.G.G., Pal, B. and Merryfield, W.J. 2011. Trends in upwelling and downwelling winds along the British Columbia shelf. Journal of Geophysical Research 116: C10023. Fraedrich, K. and Blender, R. 2003. Scaling of Atmosphere and Ocean Temperature Correlations in Observations and Climate Models. Physical Review Letters 90: 108501. Fraley, C., Leisch, F. and Maechler, M. 2006. Package ‘fracdiff’. Fractionally differenced ARIMA aka ARFIMA (p,d,q) models. R package version 1.4-2, 13 pp. Franzke, Ch.L.E. 2012. Nonlinear trends, long-range dependence, and climate noise properties of surface temperature. Journal of Climate 25: 4172-4183. 161 Franzke, Ch.L.E. 2013. Persistent regimes and extreme events of the North Atlantic atmospheric circulation. Philosophical Transactions of the Royal Society A 371: 20110471. Franzke, Ch.L.E., Graves, T., Watkins, N.W., Gramacy, R.B. and Hughes, C. 2012. Robustness of estimators of long-range dependence and self-similarity under non-Gaussianity. Philosophical Transactions of the Royal Society A 370: 1250-1267. Freeland, H.J. 1990. Sea surface temperature along the coast of British Columbia: regional evidence for a warming trend. Canadian Journal of Fisheries and Aquatic Sciences 47: 346-350. Frontier, S. 1987. Applications of Fractal Theory to Ecology. In: Legendre, P., Legendre, L., Developments in Numerical Ecology. Springer-Verlag, Berlin, pp. 335-378. Gaines, S.D. and Costello, Ch. 2013. Forecasting fisheries collapse. Proceedings of the National Academy of Sciences 110: 15859-15860. Gao, J., Hu, J., Tung, W.-W., Cao, Y., Sarshar, N. and Roychowdhury, V.P. 2006. Assessment of long-range correlation in time series: How to avoid pitfalls. Physical Review E 73: 016117. Garcia, J., Arteche, J., and Murillas, A. 2010. Fractional integration analysis and its implications on profitability: The case of the mackerel market in the Basque Country. Fisheries Research 106: 420-429. 162 Gençay, R., Selçuk, F. and Whitcher, B. 2002. An Introduction to Wavelets and other Filtering Methods in Finance and Economics, Academic Press, San Diego, California, 359 pp. Ghil, M., Yiou, P., Hallegatte, S., Malamud, B.D., Naveau, P., Soloviev, A., Friederichs, P., Keilis-Borok, V., Kondrashov, D., Kossobokov, B., Mestre, O., Nicolis, C., Rust, H.W., Shebalin, P., Vrac, M., Witt, A. and Zaliapin, I. 2011. Extreme events: dynamics, statistics and prediction. Nonlinear Processes in Geophysics 18: 295-350. Glaser, S.M., Fogarty, M.J., Liu, H., Altman, I., Hsieh, Ch.-H., Kaufman, L., MacCall, A.D., Rosenberg, A.A., Ye, H. and Sugihara, G. 2014. Complex dynamics may limit prediction in marine fisheries. Fish and Fisheries 15: 616-633. Goldberger, A. 1996. Non-linear dynamics for clinicians: chaos theory, fractals, and complexity at the bedside. Lancet 347: 1312-1314. Goldberger, A.L., Amaral, L.A.N., Hausdorff, J.M., Ivanov, P.Ch., Peng, C.-K. and Stanley, H.E. 2002. Fractal dynamics in physiology: Alterations with disease and aging. Proceedings of the National Academy of Sciences 99 (Suppl. 1): 2466-2472. Govindan, R.B., Wilson, J.D., Preibl, H., Eswaran, H., Campbell, J.Q. and Lowery, C.L. 2007. Detrended fluctuation analysis of short datasets: An application to fetal cardiac data. Physica D 226: 23-31. Gower, J.F. 2002. Temperature, wind and wave climatologies, and trends from marine metereological buoys in the Northeast Pacific. Journal of Climate 15: 3709-3718. 163 Guichard, F., Halpin, P.M., Allison, G.W., Lubchenco, J. and Menge, B.A. 2003. Mussel Disturbance Dynamics: Signatures of Oceanographic Forcing from Local Interactions. The American Naturalist 161: 889-904. Guttal, V. and Jayaprakash, C. 2008. Changing skewness: an early warning signal of regime shifts in ecosystems. Ecology Letters 11: 450-460. Hall, A. and Manabe, S. 1997. Can local linear stochastic theory explain sea surface temperature and salinity variability? Climate Dynamics 13: 167-180. Halley, J.M and Inchausti, P. 2004. The increasing importance of 1/f noises as models of ecological variability. Fluctuations and Noise Letters 4: R1–R26. Halley, J.M and Stergiou, K.I. 2005. The implications of increasing variability of fish landings. Fish and Fisheries 6: 266-276. Halley, J.M., Hartley, S., Kallimanis, A.S., Kunin, W.E., Lennon, J.J. and Sgardelis, S.P. 2004. Uses and abuses of fractal methodology in ecology. Ecology Letters 7: 254-271. Hammer, C., von Dorrien, Ch., Hopkins, Ch. C. E., Köster, F.W., Nilssen, E.M., St. John, M., and Wilson, D.C. 2010. Framework of stock-recovery strategies: analyses of factors affecting success and failure. ICES Journal of Marine Science 67: 1849-1855. Hannah, R.W. 2011. Variation in the distribution of ocean shrimp (Pandalus jordani) recruits: links with coastal upwelling and climate change. Fisheries Oceanography 20: 305-313. 164 Hanson, P.J., Vaughan, D.S., and Narayan, S. 2006. Forecasting annual harvests of Atlantic and Gulf menhaden. North American Journal of Fisheries Management 26: 753-764. Harte, D. 2001. Multifractals: theory and applications. Chapman & Hall/CRC Press, Boca Raton, Florida, 248 pp. Haslett, J. and Raftery, A.E. 1989. Space-time modelling with long-memory dependence: assessing Ireland’s wind power resource (with discussion). Applied Statistics 38: 1-50. Hastings, H.M. and Sugihara, G. 1993. Fractals: A User's Guide for the Natural Sciences. Oxford University Press, Oxford, 248 pp. Henson, S.A. and Thomas, A.C. 2007. Phytoplankton scales of variability in the California Current System: 2. Latitudinal variability. Journal of Geophysical Research 112 (C07018): 1-11. Hickey, B.M. 1998. Coastal oceanography of western North America from the tip of Baja California to Vancouver Island. In: Robinson, A.R. and K.H. Brink (Eds). The Sea, pp. 345-391. Hilborn, R. and Mangel, M. 1997. The ecological detective. Confronting models with data. Princeton University Press, Princeton, N.J. Hollister, H.J. and Sandnes, A.M. 1972. Sea surface temperature and salinity at shore stations on the British Columbia coast, 1914-1970. Fisheries and Oceans, Marine Science Directorate, Pacific Marine Sciences Report No. 72-13, Victoria, B.C. 165 Hollowed, A.B., Hare, S.R. and Wooster, W.S. 2001. Pacific Basin climate variability and patterns of Northeast Pacific marine fish production. Progress in Oceanography 49: 257-282. Hourston, R. and Thomson, R. 2013. State of physical, biological, and selected fishery resources of Pacific Canadian marine ecosystems in 2012. In Irvine, J.R. and W.R. Crawford (Eds). Physical oceanographic conditions on the shelf/shelf-break off Vancouver Island. Canadian Science Advisory Secretariat (CSAS) Research Document 2013/032, Pacific Region, Fisheries and Oceans, Canada, pp. 40-45. Hsieh, C.-H., Reiss, Ch.S., Hunter, J.R., Beddington, J.R., May, R.M. and Sugihara, G. 2006. Fishing elevates variability in the abundance of exploited species. Nature 443: 859-862. Hsieh, W.W., Ware, D.M. and Thomson, R.E. 1995. Wind-induced upwelling along the west coast of North America, 1899-1988. Canadian Journal of Fisheries and Aquatic Sciences 52: 325-334. Hu, K., Ivanov, P.Ch., Chen, Z., Carpena, P. and Stanley, H.E. 2001. Effect of trends on detrended fluctuation analysis. Physical Review E 64: 011114. Hutchings, J.A. 2000. Collapse and recovery of fishes. Nature 406: 882-885. Hutchings, J.A. and Reynolds, J.D. 2004. Marine fish population collapses: consequences for recovery and extinction risk. BioScience 54: 297-309. Inchausti, P. and Halley, J. 2001. Investigating long-term ecological variability using the Global Population Dynamics Database. Science 293: 655-657. 166 Inchausti, P. and Halley, J. 2002. The long-term temporal variability and spectral colour of animal populations. Evolutionary Ecological Research 4: 1033-1048. Ivanov, P.Ch. 2003. Long-Range Dependence in Heartbeat Dynamics. In: Processes with Long-Range Correlations. Theory and Applications. Rangarajan, G. and M. Ding (Eds.) Springer-Verlag, Berlin, pp. 339-372. Janicki, A. and Weron, A. 1994. Simulation and Chaotic Behavior of Alpha-stable Stochastic Processes. Marcel Dekker, New York, 355 pp. Jánosi, I.M. and Müller, R. 2005. Empirical mode decomposition and correlation properties of long daily ozone records. Physical Review E 71: 056126. Journel, A.G. and Huijbregts, C.J. 1978. Mining Geostatistics. Academic Press, London. Kantelhardt, J.W. 2009. Encyclopedia of Complexity and Systems Science. In: Meyers, R.A. (Ed.), SpringerScience+Business Media, LCC, New York, pp. 3754-3778. Kantelhardt, J.W. 2011. Fractal and Multifractal Time Series. In R.A. Meyers (Ed). Mathematics of Complexity and Dynamical Systems, Springer Science+Business Media, New York, pp. 463-487. Kantelhardt, J.W., Koscielny-Bunde, E., Rybski, D., Braun, P., Bunde, A. and Havlin, S. 2006. Long-term persistence and multifractality of precipitation and river runoff records. Journal of Geophysical Research 111: D01106. Kantelhardt, J.W., Koscielny-Bunde, E., Rego, H.H., Havlin, S. and Bunde, A. 2001. Detecting long-range correlations with detrended fluctuation analysis. Physica A 295: 441-454. 167 Kantelhardt, J.W., Zschiegner, S.A., Koscielny-Bunde, E., Havlin, S., Bunde, A. and Stanley, H.E. 2002. Multifractal detrended fluctuation analysis of nonstationary time series. Physica A 316: 87-114. Katul, G., Lai, Ch.-T., Schäfer, K., Vidakovic, V., Albertson, J., Ellsworth, D. and Oren, R. 2001. Multiscale analysis of vegetation surface fluxes: from seconds to years. Advances in Water Resources 24: 1119-1132. Kettani, H. and Gubner, J.A. 2006. A novel approach to the estimation of the long-range dependence parameter. IEEE Transactions on Circuits and Systems-II: Express Briefs 53: 463-467. Khaliq, M.N. and Gachon, P. 2010. Pacific decadal oscillation climate variability and temporal pattern of winter flows in Northwestern North America. Journal of Hydrometeorology 11: 917-933. Kiss, P., Müller, R. and Jánosi, I.M. 2007. Long-range correlations of extrapolar total ozone are determined by the global atmospheric circulation. Nonlinear Processes in Geophysics 14: 435-442. Kiyani, K.H., Chapman, S.C. and Watkins, N.W. 2009. Pseudononstationarity in the scaling exponents of finite-interval time series. Physical Review E 79: 036109. Kleinen, T., Held, H. and Petschel-Held, G. 2003. The potential role of spectral properties in detecting thresholds in the Earth system: application to the thermohaline circulation. Ocean Dynamics 53: 53-63. 168 Kolmogorov, A.N. 1962. A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high Reynolds number. Journal of Fluid Mechanics 13: 82-85. Koscielny-Bunde, E., Bunde, A., Havlin, S., Roman, H.E., Goldreich, Y. and Schellnhuber, H.J. 1998. Indication of a Universal Persistence Law Governing Atmospheric Variability. Physical Review Letters 81: 729-732. Koscielny-Bunde, E., Kantelhardt, J.W., Braun, P., Bunde, A. and Havlin, S. 2006. Long-term persistence and multifractality of river runoff records: Detrended fluctuation studies. Journal of Hydrology 322: 120-137. Kravchenko, A.N., Boast, Ch.W. and Bullock, D.G. 1999. Multifractal Analysis of Soil Spatial Variability. Agronomy Journal 91: 1033-1041. Kurnaz, M.L. 2004. Detrended fluctuation analysis as a statistical tool to monitor the climate. Journal of Statistical Mechanics: Theory and Experiment P07009. Lavallée, D., Lovejoy, S., Schertzer, D. and Ladoy, P. 1993. Nonlinear Variability of Landscape Topography: Multifractal Analysis and Simulation. In: Lam, N. and L. De Cola (Eds). Fractals in Geography. Prentice Hall, New Jersey, pp. 158-192. Legatt, R., Polyakov, I.V., Bhatt, U.S., Zhang, X. and Bekryaev, R.V. 2012. North Atlantic variability driven by stochastic forcing in a simple model. Tellus A 64: 18695. Lenton, T.M., Livina, V.N., Dakos, V., van Nes, E.H. and Scheffer, M. 2012. Early warning of climate tipping points from critical slowing down: comparing methods 169 to improve robustness. Philosophical Transactions of the Royal Society A 370: 1185-1204. Lewis, P.A.W. and Ray, B.K. 1997. Modeling long-range dependence, nonlinearity, and periodic phenomena in sea surface temperatures using TSMARS. Journal of the American Statistical Association 92: 881-893. Lewis, G.M., Austin, P.H. and Szczodrak, M. 2004. Spatial statistics of marine boundary layer clouds. Journal of Geophysical Research 109: D04104. Li, B.L. 2000. Fractal geometry applications in description and analysis of patch patterns and patch dynamics. Ecological Modelling 132: 33-50. Lindegren, M., Checkley Jr., D.M., Rouyer, T., MacCall, A.D. and Stenseth, N. Chr. 2013. Climate, fishing, and fluctuations of sardine and anchovy in the California Current. Proceedings of The National Academy of Sciences 110: 13672-13677. Litzow, M.A., Mueter, F. and Urban, J.D. 2013. Rising catch variability preceded historical fisheries collapses in Alaska. Ecological Applications 23: 1475-1487. Litzow, M.A. and Mueter, F.J. 2014. Assessing the ecological importance of climate regime shifts: An approach from the North Pacific Ocean. Progress in Oceanography 120: 110-119. López, J.L. and Contreras, J.G. 2013. Performance of multifractal detrended fluctuation analysis on short time series. Physical Review E 87: 022918. 170 Lovejoy, S., Currie, W.J.S., Tessier, Y., Claereboudt, M.R., Bourget, E., Roff, J.C. and Schertzer, D. 2001. Universal multifractals and ocean patchiness: phytoplankton, physical fields and coastal heterogeneity. Journal of Plankton Research 23: 117-141. Lovejoy, S. and Schertzer, D. 2012. Haar wavelets, fluctuations and structure functions: convenient choices for geophysics. Nonlinear Processes in Geophysics 19: 513-527. Lovejoy, S. and Schertzer, D. 2013. The weather and climate: emergent laws and multifractal cascades. Cambridge University Press, 475 pp. Ma, Q.D.Y., Bartsch, R.P., Bernaola-Galván, P., Yoneyama, M. and Ivanov, P. Ch. 2010. Effect of extreme data loss on long-range correlated and anticorrelated signals quantified by detrended fluctuation analysis. Physical Review E 81: 031101. MacCall, A.D. 2011. The sardine-anchovy puzzle. In: Shifting Baselines. The Past and the Future of Ocean Fisheries. Jackson, J.B.C., K.E. Alexander, and E. Sala (eds), Island Press, 47-57. Macias, D., Landry, M.R., Gershunov, A., Miller, A.J. and Franks, P.J.S. 2012. Climatic control of upwelling variability along the western North-American coast. PLoS ONE 7: e30436. MacIntosh, A.J., Pelletier, L., Chiaradia, A., Kato, A. and Ropert-Coudert, Y. 2013. Temporal fractals in seabird foraging behaviour: diving through the scales of time. Scientific Reports 3: 1884. Mahdi, E., Xiao, K.-J. and McLeod, A.I. 2014. Package ‘portes’. Portmanteau tests for univariate and multivariate time series models. R package version 2.1-3, 45 pp. 171 Makris, N.C., Ratilal, P., Symonds, D.T., Jagannathan, S., Lee, S. and Nero, R.W. 2006. Fish Population and Behavior Revealed by Instantaneous Continental Shelf-Scale Imaging. Science 311: 660-663. Mallat, S. 2009. A Wavelet Tour of Signal Processing. The Sparse Way. Academic Press, Third Edition, 805 pp. Mandelbrot, B. and Wallis, J. 1968. Noah, Joseph, and operational hydrology. Water Resources Research 4: 909-918. Mandelbrot, B.B. and Wallis, J.R. 1969. Computer experiments with fractional Gaussian noises, Parts 1, 2, 3. Water Resources Research 5: 228-267. Mandelbrot, B.B., and van Ness, J.W. 1968. Fractional Brownian motions, fractional noises and applications. SIAM Review 10: 422-437. Mandelbrot, B.B. 1977. Fractals: Form, Chance, and Dimension. Freeman W.H., & Co, San Francisco. Mandelbrot, B.B. 1983. The Fractal Geometry of Nature. Freeman W.H., & Co, San Francisco. Mandelbrot, B.B. 1985. Self-affine fractals and fractal dimension. Physica Scripta 32: 257-260. Mandelbrot, B.B. 1998. Is Nature Fractal. Science 279: 783-785. Maraun, D., Rust, H.W. and Timmer, J. 2004. Tempting long-memory - on the interpretation of DFA results. Nonlinear Processes in Geophysics 11: 495-503. 172 Marković, D. and Koch, M. 2005. Sensitivity of Hurst parameter estimation to periodic signals in time series and filtering approaches. Geophysical Research Letters 32: L17401. Martell, S.J.D. 2002. Variation in pink shrimp populations off the west coast of Vancouver Island: oceanographic and trophic interactions. Ph.D. thesis, Department of Zoology, The University of British Columbia, Vancouver, B.C. Martell, S., Boutillier, J., Nguyen, H. and Walters, C. 2000. Reconstruction of the offshore Pandalus jordani trawl fishery off the west coast of Vancouver Island and simulating alternative management policies. CSAS Research Document 2000/149, 38 pp. McClatchie, S., Greene, Ch.H., Macaulay, M.C. and Sturley, D.R.M. 1994. Spatial and temporal variability of Antarctic krill: implications for stock assessment. ICES Journal of Marine Science 51: 11-18. McCulloch, J.H. 1986. Simple consistent estimator of stable distribution parameters. Communications in Statistics - Simulation and Computation 15: 1109-1136. Ménard, F., Marsac, F. Bellier, E. and Cazelles, B. 2007. Climatic oscillations and tuna catch rates in the Indian Ocean: a wavelet approach to time series analysis. Fisheries Oceanography 16: 95-104. Mendelssohn, R. 1989. A reanalysis of recruitment estimates of the Peruvian anchoveta in relationship to other population parameters and the surrounding environment. In Pauly, D., Muck, P., Mendo, J. and I. Tsukayama (eds). The Peruvian upwelling 173 ecosystem: dynamics and interactions. ICLARM Conference Proceedings 18, Instituto del Mar del Perú (IMARPE), Callao, Perú; Deutsche Gesellschaft für Technische Zusammenarbeit (GTZ) GmbH, Eschbom, Federal Republic of Germany and International Center for Living Aquatic Resources Management (ICLARM), Manila, Philippines, pp. 364-385. Mendelssohn, R. 1990. Relating fisheries to the environment in the Gulf of Guinea: information, causality and long-term memory. In: Cury, P. and Cl. Roy (Eds). The Impact of Climatic Variations on Pelagic Fish Stocks of West Africa. ORSTOM, Dakar, pp. 446-465. Mendes, H.C., Murta, A. and Mendes, R.V. 2014. Long range dependence and the dynamics of exploited fish populations. arXiv:1406.1028v1. (Accessed 4 June 2014). Meneveau, C. and Sreenivasan, K.R. 1987. Simple multifractal cascade models for fully developed turbulence. Physical Review Letters 59: 1424-1427. Möllmann, Ch. and Diekmann, R. 2012. Marine ecosystem regime shifts induced by climate and overfishing: a review for the northern hemisphere. Advances in Ecological Research 47: 303-347. Monetti, R.A., Havlin, S. and Bunde, A. 2003. Long-term persistence in the sea surface temperature fluctuations. Physica A 320: 581-589. Monin, A.S. and Yaglom, A.M. 1975. Statistical Fluid Mechanics: Mechanics of Turbulence. MIT Press, Cambridge, Massachusetts, 874 pp. 174 Montanari, A., Taqqu, M.S. and Teverovsky, V. 1999. Estimating long-range dependence in the presence of periodicity: an empirical study. Mathematical and Computer Modelling 29: 217-228. Montes, R.M. 2004. Fractality and multifractality in the temporal dynamics of a fishery: the case of the southern hake (Merluccius australis). M.Sc. thesis, Department of Oceanography, University of Concepción, Concepción, Chile. Montes, R.M., Perry, R.I., Pakhomov, E.A., Edwards, A.M. and Boutillier, J.A. 2012. Multifractal patterns in the daily catch time series of smooth pink shrimp (Pandalus jordani) from the west coast of Vancouver Island, Canada. Canadian Journal of Fisheries and Aquatic Sciences 69: 398-413. Mudelsee, M. 2010. Climate Time Series Analyses. Classical Statistical and Bootstrap Methods. Springer Science+Business Media B.V., 474 pp. Murray, L.G., Hinz, H., Hold, N. and Kaiser, M.J. 2013. The effectiveness of using CPUE data derived from Vessel Monitoring Systems and fisheries logbooks to estimate scallop biomass. ICES Journal of Marine Science 70: 1330-1340. Myers, R.A., Hutchings, J.A. and Barrowman, N.J. 1997. Why do fish stocks collapse? The example of cod in Atlantic Canada. Ecological Applications 7: 91-106. Neubauer, P., Jensen, O.P., Hutchings, J.A and Baum, J.K. 2013. Resilience and recovery of overexploited marine populations. Science 340: 347-349. 175 Nieves, V., Llebot, C., Turiel, A., Solé, J., García-Ladona, E. and Estrada, M. 2007. Common turbulent signature in sea surface temperature and chlorophyll maps. Geophysical Research Letters 34: L23602. Niwa, H.-S. 2006. Recruitment variability in exploited aquatic populations. Aquatic Living Resources 19: 195-206. Niwa, H.-S. 2007. Random-walk dynamics of exploited fish populations. ICES Journal of Marine Science 64: 496-502. Okey, T.A., Alidina, H.M., Lo, V. and Jessen, S. 2014. Effects of climate change on Canada’s Pacific marine ecosystems: a summary of scientific knowledge. Review of Fish Biology and Fisheries 24: 519-559. Osborne, A.R. and Caponio, R. 1990. Fractal trajectories and anomalous diffusion for chaotic pattern motions in 2D turbulence. Physical Review Letters 64: 1733-1736. Osborne, A.R., Kirwan, A.D., Provenzale, A. and Bergamasco, L. 1989. Fractal drifter trajectories in the Kuroshio extension. Tellus 41A: 416-435. Oswiecimka, P., Kwapien, J., and Drozdz, S. 2006. Wavelet versus detrended fluctuation analysis of multifractal structures. Physical Review E 74: 016103. Overland, J.E., Percival, D.B. and Mofjeld, H.O. 2006. Regime shifts and red noise in the North Pacific. Deep-Sea Research I 53: 582-588. Overland, J., Rodionov, S., Minobe, S. and Bond, N. 2008. North Pacific regime shifts: definitions, issues and recent transitions. Progress in Oceanography 77: 92-102. 176 Pascual, M., Ascioti, F.A. and Caswell, H. 1995. Intermittency in the plankton: a multifractal analysis of zooplankton biomass variability. Journal of Plankton Research 17: 1209-1232. Patterson, R.M., Chang, A.S., Prokoph, A., Roe, H.M. and Swindles, G.T. 2013. Influence of the Pacific Decadal Oscillation, El Nino-Southern Oscillation and solar forcing on climate and primary productivity changes in the northeast Pacific. Quaternary International 310: 124-139. Pelletier, J.D. and Turcotte, D.L. 1999. Self-Affine time Series: II. Applications and Models. In: Dmowska, R. and B. Saltzman (eds). Advances in Geophysics. Long-range Persistence in Geophysical Time Series. Academic Press, San Diego, pp. 91-175. Peng, C.-K., Buldyrev, S.V., Goldberger, A.L., Havlin, S., Sciortino, F., Simons, M. and Stanley, H.E. 1992. Long-range correlations in nucleotide sequences. Nature 356: 168-170. Peng, C-K., Buldyrev, S., Havlin, S., Simons, M., Stanley, H. and Goldberger, A. 1994. Mosaic organization of DNA nucleotids. Physical Review E 49: 1685-1689. Peng, C.-K., Havlin, S., Stanley, H.E. and Goldberger, A.L. 1995. Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. Chaos 5: 82-87. Percival, D.B. and Walden, A.T. 2000. Wavelets Methods for Time Series Analysis. Cambridge University Press, New York, 594 pp. 177 Percival, D.B., Overland, J.E. and Mofjeld, H.O. 2001. Interpretation of North Pacific variability as a short-and long-memory process. Journal of Climate 14: 4545-4559. Perreti, Ch.T. and Munch, S.B. 2012. Regime shifts fail under noise levels commonly observed in ecological systems. Ecological Applications 22: 1772-1779. Perry, R.I., Boutillier, J.A. and Foreman, M.G.G. 2000. Environmental influences on the availability of smooth pink shrimp, Pandalus jordani, to commercial fishing gear off Vancouver Island, Canada. Fisheries Oceanography 9: 50-61. Perry, R.I., Crawford, B. and A. Sinclair. 2007. Chapter 1: Ecosystem Description. In: Lucas, B.G., Verrin, S. and R. Brown (eds). Ecosystem overview: Pacific North Coast Integrated Management Area (PNCIMA). Canadian Technical Report of Fisheries and Aquatic Sciences 2667, pp. 3-45. Petitgas, P., Secor, D.H., McQuinn, I., Huse, G. and Lo, N. 2010. Stock collapses and their recovery: mechanisms that establish and maintain life-cycle closure in space and time. ICES Journal of Marine Science 67: 1841-1848. Pilgram, B. and Kaplan, D.T. 1998. A comparison of estimators for 1/f noise. Physica D 114: 108-122. Podobnik, B., Fu, D.F., Stanley, H.E. and Ivanov, P.Ch. 2007. Power-law autocorrelated stochastic processes with long-range cross-correlations. The European Physical Journal B 56: 47-52. 178 Podobnik, B. and Stanley, H.E. 2008. Detrended Cross-Correlation Analysis: A New Method for Analyzing Two Nonstationary Time Series. Physical Review Letters 100: 084102. Pyper, B.J. and Peterman, R.M. 1998. Comparison of methods to account for autocorrelation in correlation analyses of fish data. Canadian Journal of Fisheries and Aquatic Sciences 55: 2127-2140. R Development Core Team. 2005. R: A language and environment for statistical computing. R Foundation for Statistical Computing,Vienna, Austria. ISBN 3-900051-07-0. URL http://www.R-project.org. Rea, W., Reale, M., Brown, J. and Oxley, L. 2011. Long memory or shifting means in geophysical time series? Mathematics and Computers in Simulation 81: 1441-1453. Richards, G.R. 2004. A fractal forecasting model for financial time series. Journal of Forecasting 23: 587-602. Riche, O., Johannessen, S.C. and Macdonald, R.W. 2014. Why timing matters in a coastal sea: trends, variability and tipping points in the Straight of Georgia, Canada. Journal of Marine Systems 131: 34-53. Richter, A. and Dakos, V. 2015. Profit fluctuations signal eroding resilience of natural resources. Ecological Economics 117: 12-21. Roa-Ureta, R. and Arkhipkin, A.I. 2007. Short-term stock assessment of Loligo gahi at the Falkland Islands: sequential use of stochastic biomass projection and stock depletion models. ICES Journal of Marine Science 64: 3-17. 179 Robbertse, W. and Lombard, F. 2012. On maximum likelihood estimation of the long-memory parameter in fractional Gaussian noise. Journal of Statistical Computation and Simulation 84: 902-915. Rodriguez-Iturbe, I. and Rinaldo, A. 2001. Fractal River Basins. Chance and Self-Organization. Cambridge University Press, N.Y. Rothschild, B.J. 1986. Dynamics of marine fish populations. Harvard University Press, Cambridge, Massachusetts, USA. Rouyer, T. and Fromentin, J.M. 2007. Environmental noise in spawning areas: the case of Atlantic Bluefin Tuna (Thunnus thynnus). Fisheries Oceanography 16: 202-206. Rouyer, T., Fromentin, J.-M., Ménard, F., Cazelles, B., Briand, K., Pianet, R., Planque, B. and Stenseth, N.C. 2008. Complex interplays among population dynamics, environmental forcing, and exploitation in fisheries. Proceedings of the National Academy of Sciences 105: 5420-5425. Rouyer, T., Fromentin, J.-M. and Stenseth, N. Chr. 2010. Environmental noise affects the fluctuations of Atlantic large pelagic. Progress in Oceanography 86: 267-275. Rouyer, T., Fromentin, J.-M., Stenseth, N.Chr. and Cazelles, B. 2008. Analysing multiple time series and extending significance testing in wavelet analysis. Marine Ecology Progress Series 359: 11-23. Rudnick, D.L. and Davis, R.E. 2003. Red noise and regime shifts. Deep-Sea Research I 50: 691-699. 180 Rust, H.W., Mestre, O. and Venema, V.K.C. 2008. Fewer jumps, less memory: Homogenized temperature records and long memory. Journal of Geophysical Research 113: D19110. Rutherford, D.T., Barton, L.L., Gillespie, G.E. and Boutillier, J.A. 2004. Utility of historical catch data in setting reference points for the British Columbia shrimp by trawl fishery. Canadian Science Advisory Secretariat Research Document 2004/026. Rybski, D., Bunde, A., Havlin, S., Kantelhardt, J.W. and Koscielny-Bunde, E. 2011. Detrended Fluctuation Studies of Long-Term Persistence and Multifractality of Precipitation and River Runoff Records. In: Kropp, J.P. and H.-J. Schellnhuber (eds.). In Extremis, Springer-Verlag, Berlin, pp. 217-248. Samorodnitsky, G. and Taqqu, M.S. 1994. Stable Non-Gaussian Random Processes. Stochastic Models with Infinite Variance. Chapman and Hall/CRC, Boca Raton, Florida, 632 pp. Sardà, F. and Maynou, F. 1998. Assessing perceptions: Do Catalan fishermen catch more shrimp on Fridays? Fisheries Research 36: 149-157. Scheffer, M., Carpenter, S.R., Lenton, T.M., Bascompte, J., Brock, W., Dakos, V., van de Koppel, J., van de Leemput, I.A., Levin, S.A., van Nes, E.H., Pascual, M. and Vandermeer, J. 2012. Anticipating critical transitions. Science 338: 344-348. Schertzer, D. and Lovejoy, S. 1987. Physical modeling and analysis of rain and clouds by anisotropic scaling multiplicative processes. Journal of Geophysical Research 92: 9693-9714. 181 Schertzer, D., Lovejoy, S., Schmitt, F., Chigirinskaya, Y. and Marsan, D. 1997. Multifractal cascade dynamics and turbulent intermittency. Fractals 5: 427-471. Schertzer, D. and Lovejoy, S. 1989. Nonlinear Variability in Geophysics: Multifractal Simulations and Analysis. In: Pietronero, L. (Ed.). Fractals: Physical Origin and Properties. Plenum Press, New York, pp. 49-79. Schmid, P.E. 2000. Fractal properties of habitat and patch structure in benthic ecosystems. In: Fitter, A.H. and D.G. Raffaelli (Eds). Advances in Ecological Research 30, Academic Press, pp. 339-400. Schmitt, F., Vannitsem, S. and Barbosa, A. 1998. Modeling of rainfall time series using two-state renewal processes and multifractals. Journal of Geophysical Research 103: 23181-23193. Seekell, D.A., Carpenter, S.R. and Pace, M.L. 2011. Conditional heteroscedasticity as a leading indicator of ecological regime shifts. The American Naturalist 178: 442-451. Seuront, L. 2010. Fractals and Multifractals in Ecology and Aquatic Sciences. CRC Press, Boca Raton, Florida, 335 pp. Seuront, L. 2015. On uses, misuses and potential abuses of fractal analysis in zooplankton behavioral studies: a review, a critique and a few recommendations. Physica A 432: 410-434. Seuront, L. and Lagadeuc, Y. 2001. Multiscale patchiness of the calanoid copepod Temora longicornis in a turbulent coastal sea. Journal of Plankton Research 23: 1137-1145. 182 Seuront, L. and Schmitt, F.G. 2005a. Multiscaling statistical procedures for the exploration of biophysical coupling in intermittent turbulence. Part I. Theory. Deep-Sea Research II 52: 1308-1324. Seuront, L. and Schmitt, F.G. 2005b. Multiscaling statistical procedures for the exploration of biophysical coupling in intermittent turbulence. Part II. Applications. Deep-Sea Research II 52: 1325-1343. Seuront, L., Brewer, M.C., Strickler, J.R. 2004. Quantifying zooplankton swimming behavior: the question of scale. In Handbook of scaling methods in aquatic ecology. Measurement, analysis, simulation. Edited by L. Seuront and P.G. Strutton. CRC Press, Boca Raton, Florida, pp. 333-359. Seuront, L., Gentilhomme, V. and Lagadeuc, Y. 2002. Small-scale nutrient patches in tidally mixed coastal waters. Marine Ecology Progress Series 232: 29-44. Seuront, L. and Spilmont, N. 2002. Self-organized criticality in intertidal microphytobenthos patch patterns. Physica A 313: 513-539. Seuront, L., Schmitt, F., Lagadeuc, Y., Schertzer, D. and Lovejoy, S. 1996a. Multifractal analysis of phytoplankton biomass and temperature in the ocean. Geophysical Research Letters 23: 3591-3594. Seuront, L., Schmitt, F.G., Schertzer, S., Lagadeuc, Y. and Lovejoy, S. 1996b. Multifractal intermittency of Eulerian and Lagrangian turbulence of ocean temperature and plankton fields. Nonlinear Processes in Geophysics 3: 236-246. 183 Seuront, L., Schmitt, F., Lagadeuc, Y., Schertzer, D. and Lovejoy, S. 1999. Universal multifractal analysis as a tool to characterize multiscale intermittent patterns: example of phytoplankton distribution in turbulent coastal waters. Journal of Plankton Research 21: 877-922. Shao, Y.H., Gu, G.F., Jiang, Z.Q., Zhou, W.X. and Sornette, D. 2012. Comparing the performance of FA, DFA and DMA using different synthetic long-range correlated time series. Scientifics Reports 2: 835. She, Z. and Waymire, E. 1995. Quantized energy cascade and Log-Poisson statistics in fully developed turbulence. Physical Review Letters 74: 262-269. Shen, H., Zhu, Z. and Lee, T.C.M. 2007. Robust estimation of the self-similarity parameter in network traffic using wavelet transform. Signal Processing 87: 2111-2124. Sheng, H., Chen, Y.Q. and Qiu, T.S. 2011. On the robustness of Hurst estimators. IET Signal Processing 5: 209-225. Sheng, H., Chen, Y.Q. and Qiu, T.S. 2012a. Tracking performance and robustness analysis of Hurst estimators for multifractional processes. IET Signal Processing 6: 213-226. Sheng, H., Chen, Y.Q. and Qiu, T.S. 2012b. Fractional Processes and Fractional-Order Signal Processing. Techniques and Applications. Signals and Communication Technology. In: Multifractional Processes, Springer-Verlag, London, pp. 77-92. Snover, M.L. and Commito, J.A. 1998. The fractal geometry of Mytilus edulis L. spatial distribution in a soft-bottom system. Journal of Experimental Marine Biology and Ecology 223: 53-64. 184 Sornette, D. 2000. Long-Range Correlations. In: Critical Phenomena in Natural Sciences. Chaos, Fractals, Selforganization and Disorder: Concepts and Tools. Haken, H. (Ed). Springer-Verlag, Berlin, pp. 183-198. Soutar, A. and Isaacs, J.D. 1974. Abundance of pelagic fish during the 19th and 20th centuries as recorded in anaerobic sediment off California. Fishery Bulletin 72: 257-274. Steinhaus, H. 1954. Shape and area. Colloquium Mathematicum, III: 1-13. Stephens, C., Levitus, S., Antonov, J. and Boyer, T.P. 2001. On the Pacific Ocean regime shift. Geophysical Research Letters 28: 3721-3724. Stoev, S. and Taqqu, M.S. 2004. Simulation methods for linear fractional stable motion and FARIMA using the Fast Fourier Transform. Fractals 12: 95-121. Stoev, S. and Taqqu, M.S. 2005. Asymptotic self-similarity and wavelet estimation for long-range dependent fractional autoregressive integrated moving average time series with stable innovations. Journal of Time Series Analysis 26: 211-249. Stoev, S., Taqqu, M.S., Park, Ch. and Marron, J.S. 2005. On the wavelet spectrum diagnostic for Hurst parameter estimation in the analysis of internet traffic. Computer Networks 48: 423-445. Strogatz, S.H. 2001. Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering. Westview Press, 352 pp. Sugihara, G. and May, R.M. 1990. Applications of fractals in ecology. Trends in Ecology and Evolution 5: 79-86. 185 Takalo, J., Lohikoski, R. and Timonen, J. 1995. Structure function as a tool in AE and Dst time series analysis. Geophysical Research Letters 22: 635-638. Talkner, P. and Weber, R.O. 2000. Power spectrum and detrended fluctuation analysis: Application to daily temperatures. Physical Review E 62: 150-160. Taqqu, M.S. 1987. Random processes with long-range dependence and high variability. Journal of Geophysical Research 92: 9683-9686. Taqqu, M.S. and Teverovsky, V. 1998. On estimating the intensity of long-range dependence in finite and infinite variance time series. In: Adler, R.J., Feldman, R.E., and M.S. Taqqu (eds). A Practical Guide to Heavy Tails. Statistical Techniques and Applications, Birkhäuser, Boston, pp. 177-217. Taqqu, M.S., Teverovsky, V. and Willinger, W. 1995. Estimators for long-range dependence: an empirical study. Fractals 3: 785-788. Tennekoon, L., Boufadel, M.C., Lavallee, D. and Weaver, J. 2003. Multifractal anisotropic scaling of the hydraulic conductivity. Water Resources Research 39: 1193-1204. Tennekoon, L., Boufadel, M.C. and Nyquist, J.E. 2005. Multifractal characterization of airborne geophysical data at the Oak Ridge Facility. Stochastic Environmental Research and Risk Assessment 19: 227-239. Tessier, Y., Lovejoy, S. and Schertzer, D. 1993. Universal Multifractals: Theory and Observations for Rain and Clouds. Journal of Applied Meteorology 32: 223-250. Theil, H. 1972. Statistical decomposition analysis with applications in the social and administrative sciences. North-Holland Publishing Company, Amsterdam. 186 Thomson, R.E. 1981. Oceanography of the British Columbia coast. Canadian Special Publication of Fisheries and Aquatic Sciences 56, 291 pp. Tikhonov, D.A., Enderlein, J., Malchow, H. and Medvinsky, A.B. 2001. Chaos and fractals in fish school motion. Chaos, Solitons and Fractals 12: 277-288. Torre, K., Delignières, D. and Lemoine, L. 2007. Detection of long-range dependence and estimation of fractal exponents through ARFIMA modelling. British Journal of Mathematical and Statistical Psychology 60: 85-106. Torrence, Ch. and Compo, G.P. 1998. A practical guide to wavelet analysis. Bulletin of the American Meteorological Society 79: 61-78. Tsuda, A. 1995. Fractal Distribution of an Oceanic Copepod Neocalanus cristatus in the Subarctic Pacific. Journal of Oceanography 51: 261-266. Turcotte, D.L. 1997. Fractals and Chaos in Geology and Geophysics. Cambridge University Press, New York, 398 pp. van Ness, E.H. and Scheffer, M. 2007. Slow recovery from perturbations as a generic indicator of nearby catastrophic shift. American Naturalist 169: 738-747. Vinod, H.D. 2006. Maximum entropy ensembles for time series inference in economics. Journal of Asian Economics 17: 955-978. Vinod, H.D and López-de-Lacalle, J. 2010. Maximum Entropy Bootstrap for Time Series: The meboot R Package. Journal of Statistical Software 29: 1-19. 187 Voss, R.F. 1985. Random fractal forgeries. Fundamental algorithms for computer graphics, Ilkley, U.K. Vyushin, D.I. 2009. Statistical approximation of natural climate variability. A thesis submitted in conformity with the requirements for the degree of Doctor of Philosophy. Graduate Department of Physics, University of Toronto, 217 pp. Vyushin, D.I. and Kushner, P.J. 2009. Power-law and long-memory characteristics of the atmospheric general circulation. Journal of Climate 22: 2890-2904. Vyushin, D., Mayer, J. and Kushner, P. 2009. Package ‘PowerSpectrum’ documentation. R package version 0.3, 37 pp. Wakeford, R.C., Agnew, D.J., and Mees, C.C. 2009. Review of institutional arrangements and evaluation of factors associated with successful stock recovery plans. Reviews in Fisheries Science 17: 190-222. Walters, C.J. and Martell, S.J.D. 2004. Fisheries Ecology and Management. Princeton University Press, Princeton, N.J. Wang, Y-G., Lin, Y-X. and Haywood, M.D.E. 1999. A quasi-likelihood method for fractal-dimension estimation. Mathematics and Computers in Simulation 48: 429-436. Ware, D.M. 1995. A century and a half of change in the climate of the NE Pacific. Fisheries Oceanography 4: 267-277. Webster, R. and Oliver, M. 1992. Sample adequately to estimate variograms of soil properties. European Journal of Soil Science 43: 177-192. 188 Weron A., Burnecki, K., Mercik, S. and Weron, K. 2005. Complete description of all self-similar models driven by Levy stable noise. Physical Review E 71: 016113. Wilson, P.S., Tomsett, A.C. and Toumi, R. 2003. Long-memory analysis of time series with missing values. Physical Review E 68: 017103. Witt, A. and Malamud, B.D. 2013. Quantification of long-range persistence in geophysical time series: conventional and benchmark-based improvements techniques. Surveys in Geophysics 34: 541-651. Wolter, K. and Timlin, M.S. 2011. El Niño/Southern Oscillation behaviour since 1871 as diagnosed in an extended multivariate ENSO index (MEI.ext). International Journal of Climatology 31: 1074-1087. Worm, B., Barbier, E.B., Beaumont, N., Duffy, J.E., Folke, C., Halpern, B.S., Jackson, J.B.C., Lotze, H.K., Micheli, F., Palumbi, S.R., Sala, E., Selkoe, K.A., Stachowicz, J.J. and Watson, R. 2006. Impacts of biodiversity loss on ocean ecosystem services. Science 314: 787-790. Xia, X., Lazarou, G.Y. and Butler, T. 2005. Automatic Scaling Range Selection for Long-range Dependent Network Traffic. IEEE Communications Letters 9: 954-956. Yano, J.-I., Blender, R., Zhang, Ch. and Fraedrich, K. 2004. 1/f noise and pulse-like events in the tropical atmospheric surface variabilities. Quarterly Journal of the Royal Meteorological Society 130: 1697-1721. Yu, C.X., Gilmore, M., Peebles, W.A. and Rhodes, T.L. 2003. Structure function analysis of long-range correlations in plasma turbulence. Physics of Plasmas 10: 2772-2779. 189 Yuan, N., Fu, Z. and Liu, S. 2014. Extracting climate memory using fractionally integrated statistical model: a new perspective on climate prediction. Scientific Reports 4: 6577. Zar, J.H. 1999. Biostatistical Analysis. Fourth edition, Prentice Hall, Upper Saddle River, New Jersey, USA. Zhan, H. 2008. Scaling in global ocean chlorophyll fluctuations. Geophysical Research Letters 35: L01606. Zhou, W.X. 2008. Multifractal detrended cross-correlation analysis for two nonstationary signals. Physical Review E 77: 066211. Zhu, X., Fraedrich, K., Liu, Z. and Blender, R. 2010. A demonstration of long-term memory and climate predictability. Journal of Climate 23: 5021-5029. 190 Appendix A A.1 Probability distribution of filtered sea surface temperature time series. Figure A.1 a) Filtered sea surface temperature time series for the 1940-2013 period. b) Quantile-quantile plot and c) histogram of filtered time series showing a close approximation to gaussianity. Minor deviations from gaussianity were detected in b). 191 A.2 Probability distribution of filtered wind stress time series. Figure A.2 a) Filtered wind stress time series for the 1958-1998 period. b) Quantile-quantile plot and c) histogram of filtered time series. Departures from gaussianity are evident in (b). 192 A.3 Maximum likelihood estimation of Hurst coefficients for filtered SST time series. Considering that sea surface temperature filtered time series (SSTF) was classified as fractional Gaussian noise (fGn, -1<β<1), a maximum likelihood estimator (HML) for gaussian noise of short length was used in the calculation of H (Deriche and Tewfik 1993, Pilgram and Kaplan 1998, Robbertse and Lombard 2012) according to: 𝐴. 1 𝐿 𝑥,𝛽 = −𝑛2 log (𝑥![𝑅 𝛽)]!!𝑥 − 12 log( | 𝑅(β ) |) where L(x,β) is the function to be maximized, xT correspond to the transposed SSTF time series of length n=365 days, β is the estimated LRD exponent, and |R(β)| is the determinant of the covariance matrix for unit variance (Deriche and Tewfik 1993, Pilgram and Kaplan 1998). This method is robust against noise contamination and is precise, because β values estimated for individual time series are closer to the mean of β obtained from multiple simulated fGn time series when n>256 data points (Delignieres et al. 2006). After the estimation of β, Hurst coefficients were calculated as H=(β+1)/2 (Seuront 2010). It is necessary to remark that, as the covariance matrices do not depend on time series x (equation A.1), those matrices can be pre-computed for any value of β. These latter matrices, and the MATLAB code to estimate HML were obtained from Dr. Sofiane Ramdani, University of Montpellier I, France (personal communication). Time series longer than n=365 were not analyzed using this method because it is highly time consuming due to the size of covariance matrices; consequently, its use was restricted to the analysis of annual SSTF time series. HML are shown in Fig. 3.8. 193 A.4 Filtered sea surface temperature for the 1940-2013 period. Figure A.4. Filtered sea surface temperature (SSTF) time series during the 1940-2013 period. Vertical grey lines separate periods of ten years. The period June 1997- June 1998 (highlighted with a blue bar) represents one of the two strongest El Niño events detected during the last century. Negative values of temperature were detected for this latter El-Niño event because oscillations, which include El-Niño frequency, were previously removed for the elaboration of this time series. 194 A.5 Diagram showing stages and decisions adopted for the estimation of Hurst and generalized Hurst coefficients in univariate and bivariate fractal analysis. Figure A.5 a) Detrended fluctuation analysis and maximum likelihood methods used for the estimation of Hurst coefficients (HML, HDFA). b) Steps involved in the estimation of: (i) fractal Hurst coefficient (h2), (ii) generalized fractal Hurst coefficients (h2,q), (iii) bivariate fractal coefficients (ho,f), and (iv) bivariate multifractal coefficients (ho,f(q)). Coefficients are highlighted with blue circles. Oceanographic (o) and fisheries (f) time series correspond to filtered sea surface temperature/wind stress and smooth pink shrimp total daily catch anomalies, respectively. 195 A.6 Probabilities of rejecting the null hypothesis when using a long-range dependent (LRD) model to describe a short-range dependent (AR1) process. The probabilities of erroneously using a long-range dependent (LRD) model to describe a short-range autoregressive process (AR1) were calculated using Monte Carlo simulation (Percival et al. 2001, Vyushin 2009) for time series length N=500,1500,2500…27500 (Fig. A.6). Long-range and short-range dependent time series were simulated using the Hurst coefficient (H) estimated for Amphitrite Point filtered sea surface temperature time series (H=0.87, period 1940-2013) and an AR(1) coefficient equal to 0.75. Two goodness-of-fit statistics were used to assess the adequacy of using a LRD model to fit a short-range dependence process, the Ljung-Box test and the Spectral Density Function test (Percival et al. 2001). R scripts used for the estimation of described probabilities were conducted using the Power Spectrum package (Vyushin et al. 2009) written by Dr. Dmitry Vyushin, University of Toronto, Canada (www.atmosp.physics.utoronto.ca/people/vyushin/). 196 Figure A.6 Probabilities of rejecting the null hypothesis that a long-range dependent (LRD) model is adequate for describing a realization of an autoregressive process (AR1) using the Ljung-Box statistics and Spectral Density Function test. Vertical dotted lines marked time series length equals to shrimp total daily catch (N=2568), daily mean wind stress (N=14695) and daily mean sea surface temperature (N=27010) time series, respectively.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Fractal analysis of fisheries and environmental time...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Fractal analysis of fisheries and environmental time series for the development of early warning indicators Montes-Aste, Rodrigo Marco 2015
pdf
Page Metadata
Item Metadata
Title | Fractal analysis of fisheries and environmental time series for the development of early warning indicators |
Creator |
Montes-Aste, Rodrigo Marco |
Publisher | University of British Columbia |
Date Issued | 2015 |
Description | Fractal theory has been used in aquatic sciences to quantify scale-invariant relationships under the form of scaling or power laws, which represent patterns that help to understand the complex structure of exploited marine ecosystems and their populations through time. Fisheries time series are inherently variable and complex because they represent the interaction between population dynamics, oceanographic forcing’s and anthropogenic pressure (exploitation rates). Thus, to better understand the temporal dynamics of fisheries time series under the framework of fractal theory, smooth pink shrimp (Pandalus jordani) daily catch time series, daily sea surface temperature and wind stress time series from the west coast off Vancouver Island were analyzed. The ultimate goal of this thesis was the early detection of changes and critical thresholds levels over which a fishery move to a depleted or unhealthy state. This is a highly desirable objective because it gives fishery managers time to intervene in order to avoid a fishery collapse. I identified fractal patterns in smooth pink shrimp daily catches, sea surface temperature and wind stress time series. Those patterns were quantified through the Hurst coefficient (H), which integrates time series’ long-range dependence and intermittency into a single index. Inter-annual changes in Hurst coefficients from sea surface temperature appear to be related with major oceanographic anomalous conditions such regime shifts and El Niño events. Long-range dependence and intermittency were estimated for catch time series using the fractionally differenced parameter d and the Lévy index α and then combined into a single index (1/LHE indicator). Our proposed indicator, calculated using rolling windows across the history of the fishery, tracks well the vulnerable biomass for all years, which means that it is possible to have an indicator of major changes in the biomass of Pandalus jordani at the beginning and during the progress of each fishing season. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2015-10-24 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 2.5 Canada |
DOI | 10.14288/1.0220743 |
URI | http://hdl.handle.net/2429/55105 |
Degree |
Doctor of Philosophy - PhD |
Program |
Oceanography |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2015-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/2.5/ca/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2015_november_montesaste_rodrigo.pdf [ 20.43MB ]
- Metadata
- JSON: 24-1.0220743.json
- JSON-LD: 24-1.0220743-ld.json
- RDF/XML (Pretty): 24-1.0220743-rdf.xml
- RDF/JSON: 24-1.0220743-rdf.json
- Turtle: 24-1.0220743-turtle.txt
- N-Triples: 24-1.0220743-rdf-ntriples.txt
- Original Record: 24-1.0220743-source.json
- Full Text
- 24-1.0220743-fulltext.txt
- Citation
- 24-1.0220743.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0220743/manifest