Florian Mölders
email: fmoelders@diw.de.

Ulrich Volz
email: ulrich.volz@die-gdi.de.

Diw Berlin

Birgit Schmitz

Jürgen Von Hagen

Ansgar Belke

Jörg Breitung

James Reade

Trade Creation and the Status of FTAs -Empirical Evidence from East Asia 1

Keywords: JEL Classification: F13 Free trade agreements, trade, East Asia, gravity model, panel data

HAL is

Introduction

East Asia is considered to be a latecomer to Free Trade Agreements (FTAs). By the mid-1990s the region had significantly fewer trade agreement initiatives than Europe or the Americas. Nevertheless, since 2000, many initiatives were proposed, negotiated, and even implemented. 4 FTAs with East Asian participation are concentrated within the region but have also seen extensions to extra-regional trading partners throughout the world. This introduces interesting perspectives for the analysis of trade agreements. Using trade data for 14 East Asian countries,5 this paper focuses on the impact that FTAs have on trade at the different stages that East Asian economies participate in. 6 The central part of this study is an econometric analysis using panel data with the gravity model of international trade flows. To complement the discussion about anticipatory trade effects in the dimension of preimplementation status, we control for year-to-year changes in the stages of the respective FTAs.

We use suitable econometric methods to account for the potential endogeneity problem, which occurs due to the inclusion of the FTA variable and the dubious causal relationship with the volume of trade. Additionally, we distinguish between bilateral and multilateral FTAs as these forms of trade agreements may have different trade effects, due to the expectation that either bilateral or multilateral FTAs are more likely to trigger increasing bilateral trade. If this is the case, our variables in the econometric specification will reflect this difference. As we are particularly interested in the anticipatory effects of the agreements, the included variables distinguish between all (pre-) implementation stages of bilateral and multilateral FTAs.

Our results reveal at least three interesting findings: first, there are anticipatory effects preceding the actual implementation of bilateral FTAs with East Asian participation. In particular, the stage during which bilateral agreements are (close to) being negotiated exposes significant positive trade effects. Secondly, we found no significant anticipatory trade impact for multilateral trade agreements. Thirdly, the fact that we obtain somewhat contradictory results using different regression techniquesincluding pooled ordinary least squares, generalised least squares (GLS) with random effects, Tobit regression with random effects as well as the fixed effects estimations of GLS, Poisson and negative-binomialusing identical data, highlights the need for a careful selection of the most appropriate econometric methodology.

The remainder of the paper is structured as follows: Section 2 provides an overview of trade integration and trade agreements in East Asia. Section 3 then reviews the empirical literature on trade-creating effects of FTAs, with a focus on previous studies that examine anticipatory effects. In Section 4 we introduce our econometric methodology and our data, which is followed by a presentation of our results in Section 5 and conclusions in Section 6.

Background: trade integration and trade agreements in East Asia

"When it comes to East Asian regionalism, the state of play is easily summarised -it is a mess." (Baldwin 2006: 5) Since 2000, an unprecedented surge in the development of trade agreement initiatives has occurred. Together with other Asian Pacific states, East Asia recorded the fastest growing concentration in FTAs. [START_REF] Kawai | Evolving Economic Architecture in East Asia[END_REF] and others hold multiple events responsible for the increasing number of East Asian FTAs. The main thrust behind institutionalised regional economic integration in East Asia is the ever-growing economic interdependence of the countries. It is the market-driven regionalisation [START_REF] Langhammer | The Asian Way of Regional Integration: Are There Lessons from Europe? Kiel Economic Policy Papers[END_REF]) that took off in the 1980s that was responsible for the growing interdependence in East Asia, leading to factory Asia as labelled by [START_REF] Baldwin | Managing the Noodle Bowl: The Fragility of East Asian Regionalism[END_REF]. Baldwin argues that East Asian regionalism can be classified into stages. The first, closely linked to the phenomenon of factory Asia, took off when Japan was in need of low-skilled labour due to an extraordinary increase in domestic incomes and the subsequent loss of comparative advantage in the production of labourintensive goods. This resulted in outsourcing and induced a build-up of vertical production networks throughout East Asia, accompanied by the hollowing-out of the Japanese economy.

Throughout the outsourcing period, from the mid-1980s and during the 1990s, developing East Asia7 experienced a remarkable surge in FDI inflows that considerably shaped the regional economic architecture (ADB 2006). The rising FDI inflows contributed to an increase in intra-industry and intra-regional trade, resulting from labour-seeking and component-outsourcing FDI (Kawai 2004). 8 Both forms of FDI are aimed at making use of comparative advantages and responsible for successfully networking the region. The increasing FDI inflows in East Asia "recycled" comparative advantages by upgrading the leading economies' production and allowing less developed economies to replace labourintensive industries previously dominated by the more developed countries.

Economic integration demanded further steps towards institutionalisation as weaknesses in economic cooperation became apparent during the Asian crisis. This marked a turning point as the economic interdependence became obvious in its most negative form. The lack of cooperative mechanisms that might have avoided the contagious spread of the crisis across the region became apparent. The incapability of existing trade cooperation initiatives like the ASEAN Free Trade Area (AFTA) to cushion the crisis' effects enforced the political will in the East Asian economies for further financial and trade integration and signalled that de facto economic integration would demand de jure initiatives, thereby forming a common sense of solidarity [START_REF] Kawai | Evolving Economic Architecture in East Asia[END_REF]). 9

Furthermore, the increasing regionalisation in Europe (EU) and the Americas (NAFTA, MERCOSUR) left Asian countries out of the FTA game, resulting in a fear that East Asian goods would be discriminated against in the respective regional groupings [START_REF] Kawai | Evolving Economic Architecture in East Asia[END_REF]). 10 This is consistent with the "isolation avoidance" motive where East Asian economies recognised the advantages of being institutionally connected to its trading partners [START_REF] Dent | New Free Trade Agreements in the Asia-Pacific[END_REF]. In contrast to the trade-production motives that led to the increased regionalism in East Asia, FTAs with countries outside of the region are mostly export motivated and aimed at gaining access to extra-regional markets. Accordingly, the importance of export-led growth in East Asia is stressed by [START_REF] Dobson | Deeper Integration in East Asia: Regional Institutions and the International Economic System[END_REF], who claims that the stagnation of multilateral trade liberalisation in the 1999 Seattle WTO talks caused East Asia to pursue its own intra-and extra-regional FTAs. 11 East Asian countries' trade linkages provide a heterogeneous picture of East Asia's direction of trade and its FTA policies. As [START_REF] Baldwin | Managing the Noodle Bowl: The Fragility of East Asian Regionalism[END_REF] points out, trade integration in East Asia is characterised by overlapping trade agreements via bilateral and multilateral FTAs 8 For a comprehensive overview on the effects of regional integration agreements on FDI see [START_REF] Blomström | Regional Integration and Foreign Direct Investment[END_REF]). 9 Regional initiatives have been initiated that are supposed to cushion negative shocks that could result in overlapping effects for the region. The Chiang Mai Initiative (launched in 2000) and Asian Bond Market Initiative (launched in 2003) are the most prominent ones. 10 Baldwin's (1993) domino theory helps to comprehend East Asia's fear of discrimination: the costs of nonmembership increase when trading partners (e.g. USA, Europe) enter existing trade agreements as members' trade within an agreement is treated preferentially. For reference on this issue in the context of East Asia, see [START_REF] Dent | New Free Trade Agreements in the Asia-Pacific[END_REF]. 11 [START_REF] Medvedev | Preferential Trade Agreements and their Role in World Trade[END_REF] shows that between 2000 and 2002, only 12% of total trade in East Asia was conducted under preferential trade agreements. leading to the "Noodle Bowl" syndrome. This obscure status quo displays similarities to [START_REF] Bhagwati | US Trade Policy: The Infatuation with Free Trade Areas[END_REF] "Spaghetti Bowl" that describes overlapping regulations and different rules of origin. The danger herein is that the intended trade facilitation of FTAs backfires due to multiple trade and tariff regulations that are unclear. In the case of East Asia, there are several examples of overlapping FTAs at the bilateral and multilateral level. The Sino-Thai trade relations provide a good example: China is linked with Thailand via two implemented FTAs: the bilateral People's Republic of China-Thailand Free Trade Agreement and the multilateral ASEAN-China Comprehensive Economic Cooperation Agreement. Additionally, the proposed East Asian Free Trade Area between the ASEAN+3 countries would add a third FTA between China and Thailand.

As of September 2008, there were 110 FTAs, bilateral and multilateral, connecting East Asian economies both within and outside the region (see Table 3 in the Appendix and Figure 1). Most of these agreements were at pre-implementation stage. Five (pre-) implementation stages are distinguished: proposed, framework agreement signed, under negotiation, concluded and in effect. 12 Singapore, Thailand, China and Korea stand out with the most initiatives and in that the majority of the agreements are with partners from outside the region. The other East Asian countries stand in contrast: significantly fewer FTAs that are mostly with partners from within the region. When it comes to trade agreements in the final stage (in effect), the dominance of extra-regional FTAs of the four frontrunners is weakened.

Of the in effect FTAs for China, a slight majority were extra-regional. For Thailand it is reversed: a majority are intra-regional, as the only concluded agreements outside East Asia are with Australia and New Zealand. Singapore had, as of 2008, implemented almost half of its initiatives and Korea less than a quarter.

In order to evaluate the pace of negotiations of bilateral and multilateral agreements, we have calculated the average time period (in years) that it takes for an initiative to switch from one stage to the next: bilateral agreements take about two and a half years to reach a higher stage of implementation, whereas the multilateral ones need almost three years. 13 This gives rise to the impression that bilateral FTAs are easier to negotiate and guarantee a faster implementation.

[Figure 1 and Table 1 about here.] 12 See Table 1 for the definitions of the respective stages. 13 It is worth mentioning that there are outliers in the data on bilateral FTAs. Agreements with the United States e.g. tend to be stuck in one stage for a considerable amount of time which distorts the average number of years upwards. The United States -Indonesia FTA, which stands at the proposed stage since 1997, is an example.

As of 2008, the strengthening of intra-regional trade relations in East Asia resulted in 22 intra-regional FTAs from the proposed through implemented stages. Out of these intraregional agreements, almost two-thirds were either signed or are in effect. Extra-regional

FTAs reflect an inverted picture, consisting of a mere 26% of implemented and concluded

FTAs. The comparison in the numbers between intra-and extra-regional implemented FTAs is mirrored by the trade orientation of East Asia which is characterised by high intra-regional trade.

Figure 1, which depicts the evolution of FTAs with East Asian participation, visualises the remarkable increase at the turn of the millennium. After the Asian crisis, the 1999 breakdown of WTO trade negotiations in Seattle, and China's accession to the WTO in 2001, FTA initiatives surged and a few East Asian countries set the tone. Those disappointed by the slow pace of the liberalisation process with its main trading partners through WTO negotiations in 1999 and 2003 found an alternative: establishing bilateral FTAs.

Literature review: anticipatory trade-creating effects of FTAs

Several studies look into the possible anticipatory effects, which would alter the bilateral trade volume, of trade agreements at pre-implementation stage. [START_REF] Freund | On the Dynamics of Trade Diversion: Evidence from Four Trade Blocs[END_REF] find evidence of anticipatory investments before trade agreements become effective. They attribute their findings to adjustment effects in the business environment that are needed in advance to be prepared for the period in which the FTA is implemented. On the issue of trade related anticipation, [START_REF] Eichengreen | The Role of History in Bilateral Trade Flows[END_REF] argue that FTAs are often preceded by less official agreements that lead to an orientation towards an expected formal FTA. This reorientation may also be based on the adjustment of suppliers who redirect exports in anticipation of the future agreement. In the case of the North American Free Trade Agreement (NAFTA), [START_REF] Kose | How has NAFTA Affected the Mexican Economy? Review and Evidence[END_REF] show that there was an anticipation effect before NAFTA became effective that triggered dynamic trade and financial flow effects. 14 Croce et al. (2004) argue that anticipation depends strongly on the belief that the agreement will be implemented. Hence, smooth negotiations might increase the probability of an FTA being concluded. This would suggest that an initial anticipation effect should be visible while the FTA is in the negotiation stage.

As the duration of trade negotiations varies widely, 15 depending on what trade issues are being discussed (agriculture, industry, services), a mere look at the foregoing periods does not give satisfactory insight into how the status of FTAs actually influenced bilateral trade flows. [START_REF] Magee | New Measures of Trade Creation and Trade Diversion[END_REF] investigates anticipation and lagged effects of FTAs in a panel data analysis of 133 countries between 1980 and 1998. His conclusions strongly suggest the existence of anticipatory trade effects. More precisely, trade increases, on average, by 26% in the four years prior to trade agreement implementation.

With respect to East Asia, [START_REF] Elliott | AFTA and the Asian Crisis: Help or Hindrance to ASEAN Intra-Regional Trade?[END_REF] analyse the trade creating effects of the ASEAN FTA (AFTA) before and after its implementation with a special focus on possible anticipatory trade effects. Their results indicate that no significant trade creating effects were observed in the periods following its implementation in 1992. In contrast, trade effects appeared in the period leading to the singing of AFTA. Their estimates are based on pooled data from several time periods instead of making use of panel data and its adjunct econometric features like fixed effects estimation and the possible removal of time invariant unobservable effects. Elliott and Ikemoto though employ a variable that controls for the complementarities of the respective country pair's in their factor endowments and furthermore include variables that incorporate trade-creation and -diversion effects. We omit such variables from our analysis as our focus is strictly on the anticipatory effects of trade agreements.

In accordance with these results, [START_REF] Coulibaly | Evaluating the Trade and Welfare Effects of Developing RTAs[END_REF] estimates the anticipation effects of AFTA to be positive by pointing at increasing exports and imports five years before the agreement became effective. The author not only focuses on AFTA but also refers to other RTAs and their respective anticipatory effects. Regional groupings such as the Southern African Development Community (SADC) and the Economic Community of West African States (ECOWAS) saw member trade flows increasing in the years preceding the actual implementation. Coulibaly incorporates the anticipation of trade agreements via the inclusion of a fixed lead variable meant to represent the negotiation of an agreement.

None of the aforementioned studies on anticipatory effects distinguish directly between the different stages of implementation of the respective agreements but rather focus on foregoing time periods that are not necessarily linked to any particular (pre-) implementation stage. Our study complements previous research by taking into account the degree of implementation for the respective agreements. Our goal is to contribute to the discussion on anticipatory trade effects by including a higher level of sensitivity and going beyond a mere look at preceding time periods.

Methodology and data

The gravity model is the workhorse for the analysis of international trade relations and especially the analysis of bilateral trade flows. It is based on Newton's observation that the gravity among two objects is positively proportional to their masses and adversely to the distance between them. It was introduced to economics by Tinbergen's "Analysis of World

Trade Flows" (Tinbergen 1962: 262 ff.). Further empirical work was conducted by [START_REF] Linnemann | An Econometric Study of International Trade Flows[END_REF] but theoretical justification still lagged. [START_REF] Anderson | A Theoretical Foundation for the Gravity Equation[END_REF][START_REF] Bergstrand | The Gravity Equation in International Trade: Some Microeconomic Foundations and Empirical Evidence[END_REF], among others, introduced the needed theoretical background.16

The baseline model

The basic gravity equation states that the trade between two countries is positively related to economic masses, as reflected by the countries' GDPs, and negatively by the distance between them. Economically, this means on the one hand that the more a country pair produces (the higher the product of the two GDPs), the greater the trade between the two. On the other hand, transportation costs (which increase the price that the consumer will have to pay in the importing country) will rise with increasing distance, and hence reduce bilateral trade.

In its most reduced form the gravity model includes only the respective GDPs and distance as explanatory variables for bilateral trade flows. For our model specification, we include several other variables to the basic gravity equation, so that it takes the following form:

(1) . ) ln(

) ln( ) ln( ) ln( 7 6 5 4 3 2 1 0 ijt ijt ijt jt it jt it ij ij ij ijt I WTO Pop Pop GDP GDP Language Border D T                  
For interpretative reasons, the variables are log transformed in order to mirror the percentage increase/decrease in the dependent variable due to a change in the explanatory variables. The dependent variable ijt T reflects the total bilateral trade volume in nominal USD between countries i and j (exports plus imports) in period t. 17 The explanatory variables in Equation ( 1 Variable describing FTAs between i and j in its various stages (?)

GDPitGDPjt denounces the product of the GDPs of countries i and j. The distance between i and j is labelled as Dij. In order to properly uncover the determinants of bilateral trade flows, we include several additional explanatory variables on the ground of empirical justification in the literature. 19 Sharing the same official language as well as having a common border is expected to increase bilateral trade flows. The former via trade facilitation 20 and the latter through lower transportation costs, since a common border allows goods to be traded without using a transit country through which the goods must be transferred. A factor that reflects the self-sufficiency of the respective country pairs is added with the introduction of the population variable (see [START_REF] Frankel | Regional Trading Blocs in the World Economic System[END_REF]. Countries with a higher population are expected to 17 We concentrate on nominal trade flows instead of real flows due to data availability. There are no export or import deflators on hand for the respective countries. Only a few authors record whether they use real or nominal data. Among those who do, some follow [START_REF] Rose | Do We Really Know that the WTO Increases Trade?[END_REF] and obtain real values by deflating the nominal trade flows by the US CPI. This is criticised by [START_REF] Baldwin | Gravity for Dummies and Dummies for Gravity Equations[END_REF], who state that this method is biased, as export/import deflators are needed. Following [START_REF] Frankel | Regional Trading Blocs in the World Economic System[END_REF], we include time dummies in each regression to account for effects of global inflation and growth that are common to all country pairs in any year. 18 We use distance according to the great distance formula which makes use of the countries' most populous agglomerations. 19 For a comprehensive overview on possible variables in the gravity model and further references see [START_REF] Ghosh | Are Regional Trading Agreements Trade Creating? An Application of Extreme Bounds Analysis[END_REF]. 20 [START_REF] Melitz | Language and Foreign Trade[END_REF] recently emphasised the positive influence of a common language on the volume of bilateral trade.

be less dependent on international trade as the available natural resources (land size is approximated by the product of the countries' populations) as well as the domestic market should be positively correlated with population size. The anticipated sign of this variable thus is negative.

We also include a dummy variable indicating whether both countries are members of the World Trade Organisation (WTO). Previous research suggests that joint WTO membership should increase a country pair's trade [START_REF] Rose | Do We Really Know that the WTO Increases Trade?[END_REF]Subramanian and Wei 2007;[START_REF] Helpman | Estimating Trade Flows: Trading Partners and Trading Volumes[END_REF][START_REF] Balding | Joining the World Trade Organization: What is the Impact?[END_REF]. WTO membership can be interpreted as a proxy for the countries' willingness to further integrate into the world economy. 21 Empirical findings on the positive trade effects of WTO membership is relatively consistent even though there are differences in the magnitude of the effect if one compares industrialised and developing countries (Subramanian and Wei 2007). It remains without question that, for a proper empirical analysis of trade effects, the time varying variable incorporating WTO membership ought to be included whenever a panel data set is used. Our variable is a dummy that signals whether both countries are WTO members.

The term ijt I in equation ( 1) refers to the variables describing the FTAs in its various stages and will be discussed in the following section. The cornerstone of our analysis is the analysis of trade effects during the different stages of implementation for the various trade agreements involving East Asia. The econometric specification, as well as other variables included in the regressions, plays an important role when it comes to the validity and significance of the results.

Finally, ε ijt is the error term that includes both the idiosyncratic error u ijt as well as the (unobserved) country pair-specific term a ij , being constant over time. The latter error plays a key role in the fixed effects specification discussed in section 4.3, as unobserved heterogeneity in the country pairs may be cancelled out, which reduces the endogeneity problem.

Trade data for the period 1995-2007 are drawn from the IMF's Direction of Trade Statistics (2008). GDP and population data are from IMF's World Economic Outlook (2009) database. The distance, border, and language data are from CEPII (2009). For the WTO dummy, we used the information given on the homepage of the WTO. 22 We constructed a comprehensive data set covering total trade between all East Asian countries as well as trade with 78 countries outside of the region (see table 4 in the Appendix). Regarding the extra- 21 Consider the examples of China and Vietnam. Vietnam's accession to the WTO in the beginning of 2007 was expected to boost FDI [START_REF] Hau | Economic Transition in Vietnam: Doi Moi to WTO[END_REF] similar to the boost that occurred with China's 2001 accession. 22 http://www.wto.org/english/thewto_e/whatis_e/tif_e/org6_e.htm regional partners, we included every country with at least one connection to an East Asian country via an FTA (irrespective of the status this FTA is in at period t). 23 This gives a total of 15,379 observations,24 of which 15,017 are used for the estimation. The remaining trade values are missing from the data set. Out of the available observations, 13.49% are zero observations. This is of particular importance when it comes to the Poisson and Tobit regressions as these specifications allow us to include zero observations in the analysis.

Differentiating for FTAs at different stages

Most of the empirical literature on FTAs focuses on an intuitive but superficial method to include the effects of trade agreements into the regression analysis: to include a dummy variable which takes the value 1 if an agreement is implemented to capture its trade effects from the implementation period t on, and 0 otherwise. This approach has several weaknesses, one of which we try to deal with in this study: the anticipation of agreements well in advance of the actual implementation. 25In our study, we differentiate between five stages of a trade agreement: (1) proposed,

(2) framework agreement signed, (3) under negotiation, (4) signed and (5) in effect. The respective definitions follow our data for the respective FTAs with East Asian participation, the ADB Regional Integration Center (2009). This database contains all the relevant information for any trade agreement that any East Asian economy participates in. The data covers from 1995 through the end of 2007. Therefore, we are able to source information for the regressions not only from cross-sectional but also from time variant data.

We incorporate the status variables in the regression by using a very simple but effective method: whenever an FTA initiative is proposed, negotiated or concluded in period t, a respective dummy variable reflects this status. We bring in dummies for bilateral as well as multilateral trade agreements resulting in ten dummy variables for the regression analysis.

The resulting signs and significance levels of the variables will show us which stage of the implementation process had what effect on the bilateral trade volume.

Concerning the stages of multilateral FTAs, it may occur that a country pair has overlapping FTA initiatives at various levels such as a proposed agreement as well as one that is near completion. Korea, for example, has the proposed ASEAN+3 framework and the ASEAN-Korea Comprehensive Economic Cooperation Agreement completed in 2007. In such a case, the agreement at the latest stage enters the regression analysis. Therefore, ASEAN+3 is omitted in favour of the second, "higher ranked" FTA. Since this leaves out the possibility of observing additional effects coming from the preliminary ASEAN+3 agreement, we also compute another specification that includes possible scale effects of multiple, overlapping trade agreements. These results are not reported here since the significance as well as the contribution to the discussion on anticipatory trade effects was negligible.26

Addressing the possible endogeneity of the FTA variable and other complications

When the effects of FTAs on trade flows are analysed, cautious treatment of the FTA variable is a technical imperative. The problem arises through the nature of FTAs whose establishment is based on a variety of factors, political and economic. Not only do FTAs affect bilateral trade flows, but already intense trade relations may lead to its institutionalisation, indicating that market-driven integration is followed by institutions aimed at facilitating trade flows through the removal of trade barriers. Econometrically this is problematic, as the explanatory variable must be considered endogenous as it is correlated with the error term. If an explanatory variable is correlated with the error term and unobserved effects incorporated in the error term simultaneously influence the independent and dependent variables, ordinary least squares estimates are biased. In our case, it remains unclear whether trade increases through the establishment of FTAs or vice versa. Strong trade relations could well trigger an institutionalisation of economic relations; that is, market-driven precedes institutionalised trade integration. The possibility of biased estimates (without proper adjustment of the gravity model) may be expected in the case of intra-East Asian trade where FTAs have long been neglected and market driven integration preceded institutionalisation.

The literature began to pay attention to this issue starting with [START_REF] Baier | On the Endogeneity of International Trade Flows and Free Trade Agreements[END_REF], who analyse the factors influencing the establishment of FTAs. In a subsequent study, [START_REF] Magee | Endogenous Preferential Trade Agreements: An Empirical Analysis[END_REF] shows that higher trade flows increase the probability of country pairs to establishing an FTA. Even though there is a discussion about the "natural trading partners" hypothesis [START_REF] Wonnacott | Is there a Case for Free Trade Areas[END_REF], the technical perspective of estimating FTA effects using the gravity model with endogeneity controlled specifications has long been neglected.

In general, there are two ways that the endogenous variables can be accounted for: instrumental variables (IV) estimation and the application of fixed effects. With respect to IV estimation, [START_REF] Baier | Economic Determinants of Free Trade Agreements[END_REF] investigate the determinants of FTAs adding the differences in capital/labour ratios and per capita GDPs and political variables as instruments. 27 The method of eliminating the endogeneity through IV is intuitive but in practice, this is diluted due to weak instruments and the resulting biased results (see [START_REF] Baier | Do Free Trade Agreements Actually Increase Members[END_REF]. We tried specifying the model according to IV with the inclusion of the World Bank Governance Indicators (see [START_REF] Kaufmann | Governance Matters 7: Aggregate and Individual Governance Indicators 1996-2007[END_REF], using the indicators of voice and accountability, political stability, government effectiveness, regulatory quality, rule of law, and control of corruption. 28 However, these instruments were weak, as identified by an overidentification test. 29 Baier and [START_REF] Baier | Do Free Trade Agreements Actually Increase Members[END_REF] employ another method that does not need additional explanatory variables or instruments, but rather relies on eliminating unobserved effects, such as existing intense trade relations. Fixed effects incorporate any unobserved, time invariant influences and clear out possible correlations between the regressors and the error term. If a country pair has traded intensively in the past and this remains unexplained by the standard gravity variables, then this might point towards unobserved heterogeneity, which is controlled for by using fixed effects. Colonial ties that increase the probability of an FTA between a country pair are an example of unobserved heterogeneous fixed effects that we leave of out of our econometric analysis as we focus on the fixed effects results that would remove the respective coefficient of the colonial ties variable anyway. [START_REF] Baier | Do Free Trade Agreements Actually Increase Members[END_REF] as well as [START_REF] Magee | New Measures of Trade Creation and Trade Diversion[END_REF] use fixed effects and stress their appropriateness for handling the endogeneity bias. 30 27 Their results show that differences in capital/labour ratios (making use of comparative advantages), proximity of the respective trading partners and similarity of the countries' GDP (exploiting economies of scale) among others are significant determinants of FTAs. The authors correctly predicted 97.29% of trade agreements. 28 Political factors as a driving force for FTAs has been examined by [START_REF] Mansfield | Why Democracies cooperate more: Electoral Control and International Trade Agreements[END_REF] and [START_REF] Rosendorff | Do Democracies Trade More Freely[END_REF] and their findings suggest that democratic country pairs favour bilateral trade liberalisation, as do dyads of autocracies. 29 The results are not reported here but are available from the authors upon request. 30 [START_REF] Caporale | On the Bilateral Trade Effects of Free Trade Agreements between the EU-15 and the CEEC-4 Countries[END_REF] refer to the fixed effect vector decomposition method if estimation is based on a small sample.

Other issues regarding the endogeneity bias may result from the GDP variables that are only sporadically discussed in the literature. Trade enters the GDP in the national accounts but as the gravity equation only considers bilateral trade flows, these are not expected to have a large impact on a country's GDP. Therefore, we refrain from specially treating the GDP variable in our analysis.

Besides the treatment of endogeneity for the FTA variables, a number of other issues need to be taken into consideration in order to properly estimate the model and to come up with an econometric specification accommodating the characteristics of panel data. One important problem arises from the existence of zero trade flows in the gravity equation since the gravity equation uses log-transformations. As no logarithm of zero exists, [START_REF] Frankel | Regional Trading Blocs in the World Economic System[END_REF] describes three approaches that are widely applied: the researcher can omit the zero entries and estimate the equation without those observations; but these nevertheless carry important information as the explanatory variables describe non-existing trade. This leads to a sample selection bias that can be corrected for by implementing a two-step method that first estimates the probability for two countries having zero trade and then uses these predicted values for the gravity estimates (see [START_REF] Helpman | Estimating Trade Flows: Trading Partners and Trading Volumes[END_REF]. Second, zero trade flows can be replaced with arbitrary small numbers to prevent sample selection. The obvious disadvantage is the modification of the data set. As a third approach, we can add 1 to the trade volume and then transform it into logs that then lead to zero values (in logs). Again, the drawback is the arbitrary nature of this transformation. The advantage though is that the Tobit model includes all possible observations and takes the censored nature of the data (bounded from below by zero) into account.

Santos [START_REF] Silva | The Log of Gravity[END_REF] propose using a Poisson pseudo maximum likelihood estimator that allows for zero entries in the data for the dependent variable. 31 By means of this model, the dependent variable is not transformed into logarithms and the estimated equation changes into:

(2) . ) ln( ) ln( ) ln(

7 6 5 4 3 2 1 0 ijt ijt ijt jt it jt it ij ij ij ijt I WTO Pop Pop GDP GDP Language Border D T                  
The use of the Poisson model leads over to the treatment of heteroskedastic error terms. The estimator in the standard gravity model (even controlling for fixed effects) may be biased due to measurement errors and heteroskedastic error terms. According to Santos [START_REF] Silva | The Log of Gravity[END_REF], the standard gravity model is misspecified. They favour the use of the Poisson estimation that is robust in the case of heteroskedastic error terms and measurement errors in the log-transformed variables. 32 Santos [START_REF] Silva | The Log of Gravity[END_REF]

note that even if the proportionality condition ) | ( ) | ( X Y Var X Y E  33
does not hold strictly, 34 Poisson is more adequate than models that rely on the homoskedasticity assumption. In our case, we estimate the Poisson model but our data reflects overdispersion 35 , which demands a negativebinominal distribution in order to estimate the coefficients properly (see e.g. [START_REF] Allison | Fixed-Effects Negative Binomial Regression Models[END_REF][START_REF] Soloaga | Moving Forward Faster: Trade Facilitation Reform and Mexican Competitiveness. Policy Research Working Paper 3953[END_REF]. This method basically refers to a Poisson model that accounts for overdispersion 36 while still allowing for the of zero-trade flows. 37

Empirical results

To ensure an appropriate interpretation of the various variables and highlight the partly inconsistent results of the different approaches, we present the estimates of six different specifications: pooled ordinary least squares (POLS), generalised least squares with random effects (GLS RE), Tobit regression with random effects (Tobit RE) as well as the fixed effects estimations of GLS (GLS FE), Poisson (Poisson FE) and negative-binomial (NB FE).

Bilateral as well as multilateral FTA stages are included, which yields ten dummy variables.

Overlapping multilateral agreements are treated as described above, using the highest multilateral FTA status between two countries.

Table 2 presents the respective regression results, including the results of the POLS regression, which ignores the panel characteristics of the data, to illustrate the differences between the specifications. Even though there are differences in the magnitude of the respective effects, which are mainly attributable to the respective econometric specification, the major features of the standard gravity model are confirmed throughout our analysis. The coefficients of the standard variables (distance, GDP, population, common border and common language) are largely consistent with our expectations. The distance variable which 32 An empirical application of the Poisson method for gravity models can be found in [START_REF] Magee | New Measures of Trade Creation and Trade Diversion[END_REF].

33

Y denotes the dependent variable (in our case, trade flows between i and j). X refers to the explanatory variables. 34 This would violate the equidispersion condition in the Poisson model. 35 In our case, a likelihood-ratio test rejects the hypothesis that the Poisson and negative-binomial distributions are equivalent, pointing towards overdispersion. 36 Overdisperion is characterised by an increasing ratio of the variance to the conditional mean [START_REF] Wooldridge | Econometric Analysis of Cross Section and Panel Data[END_REF]. 37 For empirical analyses employing a negative-binomial distribution, we refer to [START_REF] Soloaga | Moving Forward Faster: Trade Facilitation Reform and Mexican Competitiveness. Policy Research Working Paper 3953[END_REF] and [START_REF] Fratianni | Trade Costs and Economic Development[END_REF].

measures the great circle distance between the two nations' most important agglomerations is significantly negative for all specifications, with estimates for a trade-decreasing effect for a doubling in the distance ranging from 112% (GLS RE) and 98% (POLS) to 51% for the negative-binomial fixed effects approach. Similar results are obtained for the variables indicating a common language or the product of the respective countries' GDPs. Both coefficients are significantly positive but reflect decreasing magnitude when controlling for overdispersion and inclusion of zero-trade observations. An increase of the product of the output of the country pairs by 100% increases bilateral trade by 47%, up to 122% depending on the applied method. The coefficient estimates of the population variable are largely in line with our expectations.

[Table 2 about here.]

Turning to the effect of trade agreements, throughout the various specifications, the WTO dummy is significant at the 1% level and reveals a positive trade effect of 33.64% 38 without controlling for the panel structure and any peculiarities in the POLS specification up to 63.23% in the negative-binomial regression. This suggests a strong and highly robust effect of WTO membership in the case of bilateral trade relations with East Asia.

The FTA status dummies are of special interest as these indicate the effect on bilateral trade at the individual stage compared to a situation in which two countries do not share any common FTA status. Considering bilateral FTAs, the coefficient estimates for the proposed agreement variable are positive for all the regression methods applied. Importantly, the estimate we get with negative-binomial is not only positive but also significant, indicating the existence of an anticipatory effect. As previously noted, we consider the negative-binomial approach the most appropriate due to overdispersion in the data. We also obtain robust results

for the signed framework agreement variable, which is significantly positive for all regressions except Poisson. The coefficient of the negative-binomial signals an increase in the bilateral trade volume of 39.10%. In comparison, the fixed effect equivalent of the GLS method estimates an increase of 20.92%, which is also significant. Except for the POLS regression, the results for the under negotiation stage are insignificant, which suggests that the anticipatory effects come into play during an earlier stage. We obtain mixed results for the stages of completion variables, namely the signed and in effect variables. Regarding the former, we find a weakly significant but negative sign in the negative-binomial and a 38 The 33.64% result from calculating . 3364 . 1 29 . 0  e significant negative one in the Poisson regression (at the 5% level). According to these specifications, trade is expected to decrease by 8.33% and 12.75% respectively, which is hard to make sense of economically. The coefficients for the in effect stage are highly significant and positive in the POLS, GLS and Tobit regressions but insignificant in the econometrically more adequate Poisson and negative-binomial regressions (with the coefficient being negative for Poisson and positive for negative-binomial).

The results for the multilateral trade agreement dummies are less robust than for the bilateral FTAs. The negative-binomial regressions show almost no significant coefficients with the exception of a slightly negative coefficient for the proposed stage of multilateral FTAs. Similarly, the results for the Poisson estimates are insignificant, except for positive effects of 8.33% in the framework agreement signed stage and 20.92% for implemented multilateral FTAs, which are significant at the level of 5% and 1% respectively. The GLS and Tobit specifications reveal higher significance levels in general but also vary in the signs of their respective coefficients. Let us consider the following as an example: even though the proposed stage reveals a negative effect of 13.88% in the fixed and random effects GLS regressions, 39 the consecutive stage has the opposite trait with a slightly higher significant coefficient. In the in effect stage, multilateral FTAs on average increase bilateral trade flows by 32.31% under the Tobit specification.

The results clearly show that the magnitude and sign of the estimates depend very much on the econometric specification. This is the case for the standard gravity variables as well as for the dummy variables reflecting FTA stages. Depending on our econometric specification, the results reflect anticipatory trade effects when a bilateral FTA is right before its negotiation stage. This effect seems to be robust in our various specifications and, most importantly, the negative-binomial regressions also reflect this trait. This specification arguably has the most appropriate properties as it allows for the inclusion of zero-trade flows, can be implemented in the existence of overdispersion and allows for the estimation with fixed effects. The Poisson specification, in comparison, does not address the overdispersion in our data. The differences in the coefficients and the significance levels between the negativebinomial and the Poisson method reveal a considerable discrepancy in the two models and point towards a misspecification due to overdispersion.

The Tobit estimates also control for zero-trade with the two drawbacks being, first, a slight data-manipulation (adding 1 to each trade observation for the log-transformation), and, second, missing fixed effects. Disregarding these issues, the respective results reflect even more significant trade effects under Tobit estimates for bilateral as well as multilateral trade 39 The fixed effects approach was favoured over the random effects model by a Hausman-test.

agreements. Not only does the anticipatory effect for bilateral FTAs remain highly significant but the coefficient for implemented agreements also signals increasing trade of almost 30%.

With Tobit, even the coefficients for multilateral FTAs are significantly positive at the 1% level for multilateral agreements in framework agreement singed and the in effect stages. With this specification, and controlled for year fixed effects, implemented multilateral FTAs seem to have a stronger effect on bilateral trade flows than bilateral FTAs.

Conclusions

The economies of East Asia have only since the turn of the millennium begun to pay attention to an institutionalisation of their trade relations. This has so far resulted in a vast amount of trade initiatives, on the bilateral and multilateral level. The availability of detailed data on FTAs with East Asian participation gives us the possibility to investigate the respective trade effects of different stages of implementation.

We augment the traditional gravity variables to estimate effects coming from different sorts of trade agreements and year-to-year changes in the stages of their implementation. We use different regression techniques, including pooled ordinary least squares, generalised least squares (GLS) with random effects, Tobit regression with random effects as well as the fixed effects estimations of GLS, Poisson and negative-binomial. The results depend very much on the econometric specification. We consider negative-binomial to be the most appropriate methodology because it accounts for overdispersion in the data and allows for estimation with fixed effects. The fact that we obtain partly contradictory empirical results when using different regression techniques highlights the need not only for caution when interpreting the findings, but also for prudence in the selection of the most appropriate econometric methodology.

Our results indicate that there are indeed anticipatory effects preceding the actual implementation of bilateral FTAs with East Asian participation. Particularly the stage in which bilateral trade agreements are (close to) being negotiated exposes significant positive trade effects. The results for the effect of multilateral trade agreements are less robust and, in most cases, insignificant. This suggests that the anticipation effects of bilateral agreements are larger, possibly because the realisation of bilateral agreements is considered more realistic than that of multilateral agreements, which due to the multiple parties involved are often more complicated to negotiate. This is also reflected by the afore-mentioned average time period that it takes for FTAs to reach a consecutive stage. Bilateral trade agreements correspond to a more rapid implementation process than multilateral initiatives. Another interpretation may be the actual depth of trade integration, such that bilateral FTAs provide a more comprehensive coverage of goods than multilateral agreements. The plausibility of this argument may be reconsidered in further research which takes a closer look at the heterogeneity of the included agreements.

An extension of our work would be to replicate our analysis for other country groupings, and examine whether our finding regarding the different trade creating effects of bilateral and multilateral trade agreements in pre-implementation stages remains applicable to a different regional context. This is only a matter of data availability. 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008

Figures and Tables

#

Proposed FAS Under neg. Signed In effect Source: Own illustration based on data from ADB (2009). Note: FAS= framework agreement signed, Under neg.= under negotiation; definitions of the respective statuses are given in Table 1.

Figure 1 :

 1 Figure 1: Evolution of FTAs with and within East Asia, by status; 1995-2008

  ) are denoted as follows (with expected signs in parentheses):

	o Dij	Distance between country i and j 18 (-)
	o GDPitGDPjt	Product of the countries' GDPs (+)
	o Languageij	Common official language (+)
	o Borderij	Common border (+)
	o PopitPopjt	Product of the countries' populations (-)
	o WTOij	Both trading partners are members of the WTO (+)
	o I	ijt

Table 1 :

 1 Definition of different statuses/stages of FTAs

	Stage	Definition
	Proposed	Parties are considering a free trade agreement, establishing joint
		study groups or joint task force, and conducting feasibility studies
		to determine the desirability of entering into an FTA.
	Framework agreement	Parties initially negotiate the contents of a framework agreement,
	signed	which serves as a framework for future negotiations.
	Under negotiation	Parties begin negotiations without a framework agreement.
	Signed	Parties sign the agreement after negotiations have been completed.
	In effect	When the provisions of an FTA becomes effective, e.g. when tariff
		cuts begin.

Source:

ADB (2009)

Table 2 :

 2 Regressions by status of FTAs, Total trade *, ** and * denote significance at the 1%, 5% and 10% level; Robust standard errors in parentheses (Standard errors in Tobit and Poisson were bootstrapped); POLS=Pooled ordinary least squares, GLS=Generalised least squares, RE=random effects, FE=fixed effects; 1: value denotes within R-squared

		POLS	GLS RE	GLS FE	Tobit RE	Poisson FE	NB FE
	Method	(1)	(2)	(3)	(4)	(5)	(6)
	Constant	2.60	3.79	2.90	---	---	2.33
	Ln(Distanceij)	-0.98***	-1.12***	---	-0.65***	---	-0.51***
		(0.037)	(0.100)	---	(0.106)	---	(0.094)
	Borderij	0.82***	0.85***	---	1.05***	---	-0.26
		(0.100)	(0.273)	---	(0.274)	---	(0.576)
	Languageij	0.98***	1.21***	---	0.92***	---	0.36***
		(0.065)	(0.216)	---	(0.230)	---	(0.090)
	Ln(GDPitGDPjt)	1.22***	1.09***	0.72***	0.84***	0.73***	0.47***
		(0.009)	(0.023)	(0.039)	(0.024)	(0.063)	(0.023)
	Ln(PopitPopjt)	-0.11***	0.03	-0.89***	0.03	-0.23	-0.13***
		(0.010)	(0.031)	(0.168)	(0.024)	(0.248)	(0.031)
	WTO	0.29***	0.33***	0.36***	0.35***	0.43***	0.49***
		(0.036)	(0.036)	(0.038)	(0.051)	(0.026)	(0.067)
	Bilateral -Proposedijt	0.32**	0.07*	0.04	0.16***	0.05	0.08*
		(0.127)	(0.038)	(0.033)	(0.050)	(0.041)	(0.046)
	Bilateral -FaSijt	0.88***	0.19***	0.19***	0.29***	0.13*	0.33***
		(0.206)	(0.060)	(0.048)	(0.050)	(0.071)	(0.054)
	Bilateral -UnderNegijt	0.90***	0.05	-0.03	0.09	-0.03	0.04
		(0.152)	(0.042)	(0.029)	(0.059)	(0.065)	(0.041)
	Bilateral -Signedijt	0.49*	0.03	-0.10*	0.07	-0.12**	-0.08*
		(0.285)	(0.073)	(0.051)	(0.065)	(0.047)	(0.050)
	Bilateral -InEffect.ijt	1.12***	0.21***	0.08*	0.26***	-0.13	0.06
		(0.163)	(0.062)	(0.041)	(0.088)	(0.082)	(0.096)
	Multilat. -Proposedijt	-0.09	-0.13***	-0.13***	-0.07	0.02	-0.04*
		(0.062)	(0.035)	(0.035)	(0.053)	(0.032)	(0.020)
	Multilat. -FaSijt	0.84***	0.14***	0.14***	0.17***	0.08**	0.01
		(0.095)	(0.043)	(0.042)	(0.060)	(0.041)	(0.036)
	Multilat. -UnderNegijt	0.28***	-0.09	0.08	0.02	0.03	-0.03
		(0.091)	(0.056)	(0.056)	(0.054)	(0.032)	(0.038)
	Multilat. -Singedijt	0.92***	-0.07	-0.04	0.06*	0.02	-0.03
		(0.116)	(0.053)	(0.052)	(0.033)	(0.058)	(0.034)
	Multilat. -InEffect.ijt	1.15***	0.16**	0.01	0.28***	0.19***	-0.01
		(0.083)	(0.069)	(0.062)	(0.061)	(0.068)	(0.048)
	R²	0.76	0.75	0.25 1	---	---	---
	Observations	12,991	12,991	12,991	15,017	14,211	14,211
	Year Effects	Yes	Yes	Yes	Yes	Yes	Yes
	Note: **

Table 3 :

 3 continued Status (as of September 2008)

	Free Trade Agreements	Members	Prop.	FAS	Under neg.	Sign.	In Effect

Table 3 :

 3 continued Status (as of September 2008)

	Free Trade Agreements	Members	Prop.	FAS	Under neg.	Sign.	In Effect

Table 3 :

 3 continued Status (as of September 2008)

	Free Trade Agreements	Members	Prop.	FAS	Under neg.	Sign.	In Effect

Table 3 :

 3 continued Status (as of September 2008)

	Free Trade Agreements	Members	Prop.	FAS	Under neg.	Sign.	In Effect	Intra-EA

Table 3 :

 3 continued Status (as of September 2008)

	Free Trade Agreements	Members	Prop.	FAS	Under neg.	Sign.	In Effect	Intra-EA	Extra-EA
	Trade Preferential System of the Organization of the	Bahrain, Cameroon, Guinea, Jordan, Libyan
	Islamic Conference	Arab Jamahiriya, Pakistan, Syrian Arab
		Republic, Turkey, United Arab Emirates,
		Bangladesh, Egypt, Iran, Lebanon, Malaysia,
		Senegal, Tunisia, Uganda

Throughout this paper, we will use the terms implemented (under implementation) and in effect interchangeably.

These are: Brunei Darussalam (hereafter Brunei), Cambodia, People's Republic of China (hereafter China), Hong Kong SAR (hereafter Hong Kong), Indonesia, Japan, Republic of Korea (hereafter Korea), Lao People's Democratic Republic (hereafter Laos), Malaysia, Myanmar, Philippines, Singapore, Thailand and Vietnam.

In the subsequent analysis, we use the term FTA to describe various forms of trade enhancing initiatives such as preferential trade agreements, customs unions and of course FTAs. According to[START_REF] Vicard | On Trade Creation and Regional Trade Agreements: Does Depth matter?[END_REF], the trade creation effect of different categories of trade agreements does not vary significantly.

Developing East Asia refers to the ten ASEAN members and the NIE (Newly Industrialised Economies), namely Hong Kong, Korea and Taiwan.

In another empirical analysis onNAFTA, Krueger (1999) finds that trade relations between Canada, the United States and Mexico were affected even before the FTA became effective. According to Krueger, this was not entirely attributable to the anticipation of NAFTA but also to the liberalisation process prior to 1994 (when NAFTA was implemented) that took place in Mexico.

Dent (2006: 18) provides the example of the Japan-Mexico Economic Partnership Agreement whose negotiations took 14 rounds.

A considerable number of authors contributed to the further development and justification of the gravity model. Further inputs include[START_REF] Pöyhönen | A Tentative Model for the Volume of Trade between Countries[END_REF] and[START_REF] Helpman | Imperfect Competition and International Trade: Evidence from Fourteen Industrial Countries[END_REF]. For a comprehensive summary on the theory of the gravity model, seeDeardorff (1995) or[START_REF] Feenstra | Advanced International Trade: Theory and Evidence[END_REF] who review the relevant contributions in this field.[START_REF] Anderson | Gravity with Gravitas: A Solution to the Border Puzzle[END_REF] criticise the existing setup of the gravity equation on the grounds of missing multilateral resistance terms which shall reflect a country pairs' relative position to the rest of the world economy in addition to the bilateral resistances, such as the distance between its trading partners.

In order to address the concerns raised by[START_REF] Haveman | Trade Creation and Trade Diversion. New Empirical Results[END_REF], who argue that gravity model estimates are sensitive to the sample of countries included, we have run robustness checks with different sample selections, such as adding more countries according to the total average trade volume with East Asia or by GDP, without significantly altering our results.

This number has been calculated as follows: For each of the 14 East Asian countries, we include 78 extraregional trading partners over the period of 13 years (78×13×14=14,196 observations). There are 13 intraregional trading partners for the first East Asian economy, 12 for the second, etc., resulting in 1,183 observations ((13+12+11+10+9+8+7+6+5+4+3+2+1)×13). This adds up to 15,379 observations for East Asia's intra-and extra-regional trade

relations.25 A single dummy variable prevents the differentiation between the different types of FTAs. Trade agreements may include different kinds of trade facilitation measures, or even promote FDI or monetary cooperation among the signatories. A more sophisticated analysis should incorporate and reflect trade effects coming from different types of trade agreements in order to assign the effects to certain measures.[START_REF] Magee | New Measures of Trade Creation and Trade Diversion[END_REF], for example, distinguishes between FTAs, customs unions, and preferential trade agreements.

This method directly reflected the natural ordering of the degrees of implementation; from the proposed until the concluded stage. The five stages have been ordered in an ascending manner from 1 (proposed) to 5 (in effect). For overlapping, multilateral FTAs, we summed the numbers to incorporate scale effects coming through multiple trade agreements. The results are available from the authors upon request.

For further justification of the Poisson model in non-count models, see[START_REF] Winkelmann | Econometric Analysis of Count Data[END_REF] and[START_REF] Wooldridge | Econometric Analysis of Cross Section and Panel Data[END_REF] who recommend Poisson due to its properties under the existence of heteroskedasticity.

Japan-Thailand Economic Partnership AgreementJapan, Thailand √ √ Japan-Vietnam Economic Partnership Agreement Japan, Vietnam √ √

New Zealand, Singapore √ √

China, Republic of South Africa √ √

Thailand-Peru Free Trade AgreementThailand, Peru √ √

United States-Thailand Free Trade Agreement United States, Thailand √ √ Source: Own compilation based on data from ADB (2009) Note: Prop.=Proposed; FAS=Framework agreement signed; Under. neg.=Under negotiation; Sign.=Signed; Intra-EA describes FTA initiatives that consist solely of East Asian economies. Extra-EA includes at least one country from outside of East Asia.

Appendix


Christian Guill
email: guill@fkp.tu-darmstadt.de

Barbara Drossel
email: drossel@fkp.tu-darmstadt.de

Wolfram Just
email: w.just@qmul.ac.uk

Eddy Carmack
email: eddy.carmack@dfo-mpo.gc.ca

A three-species model explaining cyclic dominance of pacific salmon

Keywords: cyclic dominance, population oscillations, single-generation cycles, Neimark Sacker bifurcation, strong resonance

The four-year oscillations of the number of spawning sockeye salmon (Oncorhynchus nerka) that return to their native stream within the Fraser River basin in Canada are a striking example of population oscillations. The period of the oscillation corresponds to the dominant generation time of these fish.

Various -not fully convincing -explanations for these oscillations have been proposed, including stochastic influences, depensatory fishing, or genetic effects. Here, we show that the oscillations can be explained as an attractor of the population dynamics, resulting from a strong resonance near a Neimark Sacker bifurcation. This explains not only the long-term persistence of these oscillations, but also reproduces correctly the empirical sequence of salmon abundance within one period of the oscillations. Furthermore, it explains the observation that these oscillations occur only in sockeye stocks originating from large oligotrophic lakes, and that they are usually not observed in salmon species that have a longer generation time.

Introduction

In many ecological systems distinct population oscillations are known, with their specific dynamical characteristics often captured by simple generic models. Among these are the spatial synchronisation of the lynx-hare oscillations in Canada [START_REF] Elton | The ten-year cycle in numbers of lynx in Canada[END_REF][START_REF] Odum | Fundamentals of Ecology[END_REF][START_REF] Blasius | Complex dynamics and phase synchronization in spatially extended ecological systems[END_REF], the chaotic oscillations of boreal rodents in Fennoscandia [START_REF] Hanski | Population oscillations of boreal rodents -Regulation by mustelid predators leads to chaos[END_REF], or the cyclic outbreak dynamics of the spruce budworm [START_REF] Ludwig | Qualitative analysis of insect outbreak systems -spruce budworm and forest[END_REF][START_REF] Royama | Population dynamics of the spruce budworm choristoneura-fumiferana[END_REF]. The four-year oscillations of sockeye salmon (Oncorhynchus nerka) in the Fraser River basin in Canada are another well-documented example of such large-scale population oscillations [START_REF] Ricker | Cycle dominance among the Fraser sockeye[END_REF][START_REF] Townsend | Population cycles in freshwater fish[END_REF][START_REF] Ricker | Cycles of abundance among Fraser River sockeye salmon (Oncorhynchus nerka)[END_REF]. They are a typical example of single-generation cycles [START_REF] Gurney | Fluctuation periodicity, generation separation, and the expression of larval competition[END_REF][START_REF] Murdoch | Single-species models for many-species food webs[END_REF]. Every fourth year (which is the dominant generation time) the abundance of these fish is at very high levels, reaching several million fish in some spawning populations, but drops to numbers between several hundred and a few ten thousand individuals in the following years (hence the term cyclic dominance). Different stocks can have their population maximum in different years (Fig. 1 andAppendix). The oscillations were reported as early as the 19th century and are evident for instance in the extremely high catches by fisheries every fourth year [START_REF] Rounsefell | The salmon and salmon fisheries of Swiftsure Bank, Puget Sound, and the Fraser River[END_REF]. This both economically and conservationally important phenomenon has been ascribed either to transient effects or to stochastic influences [START_REF] Myers | Simple dynamics underlie sockeye salmon (Oncorhynchus nerka) cycles, Can[END_REF], to depensatory predation [START_REF] Larkin | Simulation studies of the Adams River sockeye salmon, Oncorhynchus nerka[END_REF], to fishing [START_REF] Walters | Evidence against the existence of cyclic dominance in Fraser River sockeye salmon (Oncorhynchus nerka), in: Sockeye salmon (Oncorhynchus nerka) population biology and future management[END_REF], or to genetic effects [START_REF] Levy | Review of proposed mechanisms for sockeye salmon population cycles in the Fraser River[END_REF][START_REF] Walters | Genetic models for cyclic dominance in sockeye salmon (Oncorhynchus nerka), Can[END_REF], but all of these explanations are still not fully convincing [START_REF] Levy | Review of proposed mechanisms for sockeye salmon population cycles in the Fraser River[END_REF][START_REF] Ricker | Cycles of abundance among Fraser River sockeye salmon (Oncorhynchus nerka)[END_REF]. Total escapement (number of male and female adults that escape the fisheries) of the late Shuswap Lake run. Shuswap lake is one of the most important rearing lakes of juvenile sockeye salmon. Data for 5 other lakes are provided in the Appendix.

The sockeye salmon return to spawn in their native stream or lake in late summer and then die, which means that the salmon generations do not overlap. The hatched fry migrate downstream in the following spring and feed for one season in large freshwater lakes, before they migrate to the ocean, where they spend the next two and a half years. The carcasses of the adult salmon are decomposed and provide an important phosphorus input into the rearing lakes of the fry. Since the oscillations of sockeye salmon originating from different lakes are not in phase, we have clear evidence that the relevant processes causing the phenomenon of cyclic dominance take place in the rearing lakes rather than in the ocean.

Model and results of computer simulations

It is our aim to develop a generic model capturing only the essential mechanisms required for the occurrence of cyclic dominance. As such, the model is kept very simple and is not designed to quantitatively predict the population dynamics of all species in the corresponding lake ecosystems.

The model uses standard continuous population dynamics equations for sockeye fry, s n (t), their predator (e.g. rainbow trout), p n (t), and their zooplankton food (mainly daphnia), z n (t), during the growth season from spring (t = 0) to fall (t = T ) in year n, combined with a rule for calculating the three population sizes at the beginning of the next season as a function of the population sizes at the end of the previous season(s). The sockeye fry population at the end of a season gives rise to the number of spawning adults 3 years later, which in turn determines the number of sockeye fry in the following spring, and represents a nutrient input that affects the carrying capacity of the zooplankton. A small fraction of the salmon stays in the ocean for one more year and mature at age 5, thus causing a mixing between the four brood lines of a spawning population. An even smaller fraction matures at age 3, but since these fish are predominantly small-sized males [START_REF] Ricker | Cycles of abundance among Fraser River sockeye salmon (Oncorhynchus nerka)[END_REF] that do not influence the number of fertilised eggs, they are neglected in this study.

The general structure of the model, which combines continuous population dynamics during the growth season with a discrete update scheme capturing the seasonally driven migration and reproduction of the species of interest, has also been applied in a modelling study of the dynamics of fish preying on zooplankton [START_REF] Persson | Ontogenetic scaling of foraging rates and the dynamics of a size-structured consumer-resource model[END_REF]. Here, however, we have chosen a much simpler form of the continuous dynamics which is sufficient to generate cyclic dominance. It is given by the following equations of motion:

d dt s n (t) = λ • a sz z n (t) • s n (t) 1 + c s • s n (t) + z n (t) -a ps s n (t) • p n (t) 1 + c p • p n (t) + s n (t) -d s • s n (t) d dt z n (t) = z n (t) • 1 - z n (t) K n (1) -a sz z n (t) • s n (t) 1 + c s • s n (t) + z n (t) d dt p n (t) = λ • a ps s n (t) • p n (t) 1 + c p • p n (t) + s n (t) -d p • p n (t)
The variables s n , z n , and p n are biomass densities that are made dimensionless. The parameters in these equations and their numerical values used in the computer simulations are as follows. λ = 0.85 denotes the assimilation efficiency of ingested prey biomass of carnivores [START_REF] Yodzis | Body size and consumer-resource dynamics[END_REF]. The maximal per unit mass ingestion rates of salmon, a sz = 10, and predators, a ps = 1.6, scale allometrically with bodymass according to a power law with exponent -1/4 [START_REF] Brown | Toward a metabolic theory of ecology[END_REF]. We assume the predators to be approximately 1500 times heavier than the sockeye fry, which fixes the ratio a sz /a ps . The parameters d s = 1 and d p = 0.16 represent biomass loss due to respiration and mortality; they also scale allometrically with body mass.

The feeding terms include saturation at high prey densities, and a predator interference term in the denominator with interference strengths set to c s = 1 and c p = 0.2. This Beddington functional response [START_REF] Beddington | Mutual interference between parasites or predators and its effects on searching efficiency[END_REF] has been chosen since it describes consumer-resource interactions more accurately than the more frequently used Holling type II form [START_REF] Skalski | Functional responses with predator interference: viable alternatives to the Holling type II model[END_REF]. Similar predator interference terms are also used by other modellers [START_REF] Walters | Adding realism to foraging arena predictions of trophic flows in Ecosim ecosystem models: shared foraging arenas and bout feeding[END_REF].

The matching conditions used to determine the biomasses of the species at the beginning of the next season from their values at the end of the previous season(s) are given by

s n+1 (0) = γ((1 -)s n-3 (T ) + s n-4 (T )) z n+1 (0) = K n+1 (2) p n+1 (0) = p n (T ) ,
with the proportion of surviving sockeye that return to their native lakes at the age of 5 to spawn and die. The parameter γ summarises survival from smolt to adult fish (including ocean survivability and loss to fisheries), spawning success, and egg to fry survival. The various factors are estimated from empirical data [START_REF] Walters | Evidence against the existence of cyclic dominance in Fraser River sockeye salmon (Oncorhynchus nerka), in: Sockeye salmon (Oncorhynchus nerka) population biology and future management[END_REF][START_REF] Hume | Juvenile sockeye rearing capacities of three lakes in the Fraser River system[END_REF][START_REF] Pauley | Species profiles: life histories and environmental requirements of coastal fishes and invertebrates (Pacific Northwest) -Sockeye salmon[END_REF] and yield γ = 0.4.

According to empirical observations [START_REF] Levy | Review of proposed mechanisms for sockeye salmon population cycles in the Fraser River[END_REF], the zooplankton level at the end of one year has no effect on the following year. However, its carrying capacity K n+1 is a function of the nutrients provided by the adult salmon and thus a function of the initial fry biomass of that year [START_REF] Hume | Sockeye fry, smolt, and nursery lake monitoring of Quesnel and Shuswap lakes in 2004[END_REF]:

K n+1 = K 0 + κ s n+1 (0) κ 0 + s n+1 (0) . ( 3
)
(For a review on the importance of salmon-derived nutrients for freshwater ecosystems, see [START_REF] Gende | Pacific salmon in aquatic and terrestrial ecosystems[END_REF] and citations therein.) Consumer dependent resource productivity is certainly a rather uncommon phenomenon. In the lake ecosystems under consideration, it provides an important positive feedback to the system, but it does not seem to be crucial for the mathematical mechanism underlying the oscillations [START_REF] Guill | Coexistence of two periodic attractors: the degenerate Neimark Sacker bifurcation as a generic mechanism[END_REF]. The three parameters determining the carrying capacity K n of the zooplankton in year n are set to K 0 = 15, κ = 10, and κ 0 = 2. The numerical values of these parameters are difficult to determine from empirical data. However, since the nutrient concentration in the brood lakes nearly doubles in years following dominant sockeye returns compared to years following non-dominant sockeye returns [START_REF] Hume | Sockeye fry, smolt, and nursery lake monitoring of Quesnel and Shuswap lakes in 2004[END_REF], it seems reasonable that the constant fraction of the carrying capacity, K 0 , and the maximum of the variable fraction, κ, are of the same order of magnitude. By additionally increasing K 0 and c p , their order may even become reversed. This is also observed in the empirical data of some stocks.

When a parameter is varied, for instance the constant fraction of the zooplankton carrying capacity, K 0 , the dynamical pattern may change. Fig. 3 shows the biomass s n (T ) of the sockeye fry at the end of the season for different values of K 0 , from year t =5000 to year 5200. For small K 0 , all sockeye lines are equally strong (only one point is visible for each value of K 0 ). This means that the dynamics reaches a fixed point and that there is no cyclic dominance. With increasing K 0 , the fixed point eventually becomes unstable. First a bifurcation to quasiperiodic behaviour occurs (the attractor contains infinitely many points for each value of K 0 ), and then the frequency of the oscillation becomes locked at 4. Now only four points are visible for each value of the bifurcation parameter, indicating periodic oscillations.

Bifurcation analysis

In order to understand and interpret these observations, we first note that the continuous population dynamics during the season, together with the matching conditions applied between two seasons, can be viewed as a discrete map (a so-called Poincaré map), giving the sockeye and predator biomasses at the end of one year as function of the biomasses at the end of the previous years.

To obtain this map, one first has to integrate the dynamical equations over one season, giving s n (T ) and p n (T ) as a function of s n (0) and p n (0). The zooplankton can be eliminated since the initial condition of the zooplankton depends on s n (0) only. Next, one expresses s n+1 (0) and p n+1 (0) as a function of s n-3 (T ), s n-4 (T ) and p n (T ), using the matching conditions.

The mechanism which generates the population oscillation is based on the nature of the instability of the stationary state s n (T ) = s * (T ), p n (T ) = p * (T ) of our system. The corresponding bifurcation can be investigated in terms of a linear stability analysis. Close to the fixed point, the dynamics can be approximated by linear terms. Denoting the distance of the biomasses from their fixed point value by δs n = s n (T )s * (T ) etc, the linear approximation of this map has the form

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ δs n+1 δs n δs n-1 δs n-2 δs n-3 δp n+1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 0 0 m s 1 -a 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 b c m p ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ δs n δs n-1 δs n-2 δs n-3 δs n-4 δp n ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (4)
with positive parameters m s , m p , a, b, c, 1 . The general structure of the matrix in Eq. ( 4) is determined only by the matching conditions and is valid for any time-continuous model applied during the seasons.

The equations of motion of the latter only determine the numerical values of the parameters in Eq.

(4).

The first line of the matrix describes the influence on the sockeye fry of year n + 1 of the sockeye fry of year (n -3) and (n -4) (which are the parents of the fry in year (n + 1)), and of the predator population. m s and 1 are positive, since more parents imply more offspring. -a is negative, since more predators imply less fry. The other nontrivial line of this matrix, the last line, describes the influence on the predators of year (n + 1) of the sockeye fry of year (n -3) and (n -4) (which are the parents of the predator's food), and of the predator population in the previous year. All three parameters are positive, since more food implies more predator growth and since more predators in one year give rise to more predators in the next year.

The eigenvalues of this matrix determine the nature of the dynamics near the bifurcation. When all eigenvalues have an absolute value smaller than 1, the fixed point is stable, and the dynamics converges to this fixed point. When the absolute value of one or more eigenvalues is larger than 1, the fixed point is unstable, and the dynamics approaches a different attractor.

In order to understand the dynamics resulting from this matrix, we first consider the case that the matrix elements 1 , b, and c vanish. This means that all salmon return at age 4, and that the trout have a good choice of other food so that their dynamics is independent of that of the salmon fry. In this case the eigenvalues of the matrix are m p , 0, and the four fourth roots of m s . Since the four salmon lines are independent from each other in this case, the sequence δs n has trivially the period 4 and simply iterates the initial four values, with an amplitude decreasing for m s < 1 and increasing otherwise (and with the trout being independent of the salmon). When m s is increased from a value smaller than 1 to a value larger than 1, all four eigenvalues m 1/4 s cross the unit circle simultaneously, and the fixed point becomes unstable. This degeneracy is lifted when the parameters 1 , b, and c are made nonzero. As long as these parameters are not large, one can expect the four main eigenvalues to remain close to the real and imaginary axis, respectively, implying a (possibly damped) oscillation with a period close to 4.

The type of bifurcation that occurs when the fixed point becomes unstable depends on which eigenvalue crosses first the unit circle as a control parameter is increased. In order to determine the type of the bifurcation, we evaluated the parameters of the matrix numerically from our computer simulations of the time continuous model when going through the bifurcation (i. e., when K 0 increases from 6.4 to 7.5), and we calculated the eigenvalues λ of the matrix. Figure 4 shows that the bifurcation is caused by a pair of complex conjugate eigenvalues crossing the unit circle, indicating a Neimark Sacker bifurcation (the discrete version of the Hopf bifurcation). Since the dominant eigenvalues (depicted with the symbol +) are close to ±i, the period of the resulting oscillation is close to 4. In the linear, time discrete model (Eq. ( 4)), the bifurcation is mainly driven by the parameters m s and, to a minor extent, a, which correspond to reproduction of the four-year old sockeyes and predation by the rainbow trouts, respectively. The fixed point becomes unstable when either of these parameters increases. Now, it is known from the theory of bifurcations that if the period at such a bifurcation is close to 4, there occurs a strong resonance, which means that the period becomes locked exactly at the value 4 not far beyond the bifurcation point. In contrast to conventional mode locking that gives rise to Arnold tongues, a strong resonance is due to nonlinearities that are of the same order as the leading nonlinearity (which ensures the stability of the quasiperiodic orbit in the first place), and frequency locking therefore occurs over a much wider range of parameters [START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF].

The two complex conjugate eigenvalues are the first ones to cross the unit circle when a and 1 are small non-negative numbers and when m p is well below 1. Since 1 , b, and c are small compared to m s , the period of the resulting oscillation remains close to 4, and the locking at period 4 therefore occurs not far beyond the bifurcation. When these conditions are not met, the unit circle is typically first crossed at -1, and a period-doubling (or flip) bifurcation occurs.

Discussion

The four-year oscillations of sockeye salmon in the Fraser River are a typical example of singlegeneration cycles (as opposed to predator-prey cycles [START_REF] Gurney | Fluctuation periodicity, generation separation, and the expression of larval competition[END_REF]), a phenomenon that has also been extensively studied in the context of, e.g., host-parasitoid [START_REF] Godfray | Discrete and continuous insect populations in tropical environments[END_REF] or daphnia-zooplankton [START_REF] De Roos | A size dependent predator-prey interaction: who pursues whom?[END_REF] systems. In the latter case, individuals of different age or size of a consumer species directly compete with each other for a resource, while the biological situation in the present study is somewhat different: The salmon fry of different brood lines do not compete with one another directly, since they populate the rearing lakes at distinct, non-overlapping times. The interaction between the brood lines is mediated by the long-living predators, which have a net effect on the salmon that is similar to competition (apparent competition). In contrast to other models that produce single-generation cycles, the cycle period 4 observed in our system is not identical to the mean generation time, which is larger than 4 years, and not all salmon mature and reproduce at the same age. The locking of the period at 4 is due to a strong resonance, a phenomenon which plays no role in the cited models.

The three-species model analysed in this study is in principle also capable of showing conventional predator-prey cycles. If the parameters were chosen such that the complex conjugate pair of eigenvalues of the matrix in Eq. ( 4) that is close to the real axis at +1 is the first to cross the unit circle, oscillations with a period much larger than four years occur and the year-to-year variation between brood lines is rather small (results not shown). However, this type of oscillations is not seen in the sockeye populations rearing in the Fraser River basin.

The considerations in the previous section lead to two basic conditions under which strong resonance and thus period-four generation cycles can be observed: First, an increase in the number of spawning salmon must lead to a sufficiently strong increase in this number four years later (i.e. m s must be large enough). If the number of fry migrating to the ocean is dominated by other factors, such as a strongly limited carrying capacity for the fry, the population will be at a fixed point rather than on the oscillating side of the bifurcation. Second, the four salmon lines must be coupled in order for the Neimark Sacker bifurcation to occur, rather than a period doubling bifurcation. In our model, this coupling is due to a fraction of sockeye returning at age 5 instead of age 4, and, more importantly, due to the predator being sufficiently strongly coupled to sockeye dynamics.

Previous studies of salmon dynamics, used in salmon management, are based on the Ricker model [START_REF] Ricker | Cycle dominance among the Fraser sockeye[END_REF][START_REF] Hume | Juvenile sockeye rearing capacities of three lakes in the Fraser River system[END_REF] or the Larkin model [START_REF] Larkin | Simulation studies of the Adams River sockeye salmon, Oncorhynchus nerka[END_REF][START_REF] Martell | Retrospective analysis of harvest management performance for Bristol Bay and Fraser River sockeye salmon[END_REF][START_REF] Marsden | Retrospective bioeconomic analysis of Fraser River sockeye salmon fishery management[END_REF]. Because both models have a strongly limited carrying capacity for the fry, and because they do not include explicitly the predator dynamics, the only bifurcation occurring in those models is the flip bifurcation. Building on both modelling approaches, Myers et al. [START_REF] Myers | Simple dynamics underlie sockeye salmon (Oncorhynchus nerka) cycles, Can[END_REF] have demonstrated that an unstable (decaying) period-4 oscillation can be excited by stochastic driving. However, it remains unclear over how long the oscillations can be maintained in that study without a change in the phase. Furthermore, that stochastic model produces only episodically the typical sequence of a dominant, a subdominant and two weak years described in section 2 and seen in our simulations, Fig. 2).

In principle, the coupling between the four salmon lines can also occur through the food of the sockeye fry. However, there is no empirical evidence that there is a negative effect of a strong sockeye year on the daphnia population in the following year that limits the growth of the sockeye fry in that year [START_REF] Levy | Review of proposed mechanisms for sockeye salmon population cycles in the Fraser River[END_REF].

The two conditions for a strong resonance of period 4 fit together with the empirical observation that it occurs only in large oligotrophic lakes, such as those of the Fraser River basin. The smaller and ultra-oligotrophic lakes in the outer coast regions of British Columbia do not show these oscillations.

Their nutrient content is neither large enough to raise large fry numbers, nor to allow for salmon predators to become strong. On the other hand, cyclic dominance cannot be expected in nutrientrich lakes, because the spawning adults would not be an important nutrient input, and because the coupling to the predator would not be strong enough in a situation where there were more predator species and more prey species for these predators.

Since there are large fluctuations in the proportion of fry that survive to return to their nursery lake, we included noise in the matching condition for the sockeyes, Eq. ( 2), in order to determine with how much noise superimposed on the deterministic dynamics the period-4 oscillation can persist. Fig. 2b shows a data series generated with 50 percent noise in ocean survivability superimposed to the deterministic dynamics. The oscillation is still clearly visible, although the system has a phase shift every 300 years on average. With less noise, the phase shift occurs less often, and with 100 percent noise, the cyclic dominance vanishes in our simulations. When a large perturbation acts only for a limited time, the oscillation quickly recovers afterwards. In fact, the recovery following the blocking of the Fraser River migration routes early in the last century, can be seen in the non-dominant brood lines of the Shuswap stock in Fig. 1 and in the data of some of the stocks presented in the Appendix.

Unfortunately, the year 2009 has seen another large perturbation, with most Fraser sockeye expected to return that year not surviving in the ocean, so that the expected strong escapement did not occur.

The age composition of sockeye stocks was also found to strongly influence the potential of the system to show a strong resonance. The fraction of adult sockeyes that returns at the age of 5 years instead of 4 years, , was set to 0.1 in the simulations and we could show that the resonance appears for 0 < ≤ 0.3, but for very small values of (< 0.02) the non-dominant lines disappear completely with increasing K 0 . The resonance condition was best met at ≈ 0.2, where the resonance occurred nearly directly after the Neimark Sacker bifurcation. The parameter is also measured in the real populations [START_REF] Healey | The adaptive significance of age and size at maturity in female sockeye salmon (Oncorhynchus nerka)[END_REF]. For the sockeye salmon populations of the Fraser River it is approximately 0.08, while for the less productive outer coast lakes of British Columbia (where cyclic dominance is not observed) it is between 0.56 and 0.76. Some of the sockeye populations of the Bristol Bay area (Alaska), most notably the Kvichak River stocks, also exhibit strong oscillations [START_REF] Fair | Critical elements of Kvichak River sockeye salmon management, Alaska Fish[END_REF][START_REF] Rogers | Asynchrony in population dynamics of sockeye salmon in southwest Alaska[END_REF] despite a broader distribution of the age at spawning than in the Fraser River stocks [START_REF] West | Abundance, age, sex, and size statistics for pacific salmon in Bristol Bay[END_REF][START_REF] West | Abundance, age, sex, and size statistics for pacific salmon in Bristol Bay[END_REF]. However, the oscillations are not as regular, with maxima of the populations occurring every fourth or fifth year. This is consistent with our model, as it indicates a quasiperiodic oscillation rather than a fixed oscillation period associated with a strong resonance.

Our results do not rule out additional mechanisms such as depensatory fishing [START_REF] Walters | Evidence against the existence of cyclic dominance in Fraser River sockeye salmon (Oncorhynchus nerka), in: Sockeye salmon (Oncorhynchus nerka) population biology and future management[END_REF] or genetic effects [START_REF] Walters | Genetic models for cyclic dominance in sockeye salmon (Oncorhynchus nerka), Can[END_REF], which would reduce the population sizes of weak lines to values smaller than those resulting from our model. However, these additional assumptions are not needed to explain the occurrence of cyclic dominance in the first place.

Appendix: Escapement data for sockeye stocks showing cyclic dominance

Many spawning stocks of sockeye salmon in the upper Fraser River basin are not stationary but have been dramatically increasing in size over the last decades (e.g. Figure A1 a, b, and d), following a massive disturbance of the stocks at the beginning of the 20th century. Nevertheless, the oscillatory pattern has clearly emerged again.

In some stocks (most notably the Bowron Lake stock, Figure A1 c), cyclic dominance is observed only episodically. This dynamical pattern may appear when the system is close to the bifurcation point.

The sockeye salmon stock of Seton Lake (Figure A1 e) also displayed cyclic dominance, but for nearly two decades (five cycles), the dominant line is followed by the weakest line instead of the sub-dominant line. In this stock, cyclic domaince broke down in the early 1990's.

Figure 1 :

 1 Figure1: Empirical time series of sockeye salmon abundance. Total escapement (number of male and female adults that escape the fisheries) of the late Shuswap Lake run. Shuswap lake is one of the most important rearing lakes of juvenile sockeye salmon. Data for 5 other lakes are provided in the Appendix.

Fig. 2a shows

  Fig.2ashows a time series of the biomass of the sockeye fry at the end of the growth season. The first 300 years are cut off to show only the stable periodic oscillation with one strong year followed by one intermediate year and two weak years, just as in the empirical data of most sockeye stocks showing cyclic dominance (Fig.1 and Appendix). When the simulation parameters are chosen differently (e.g., c s = 0.2), the difference between the strong and the intermediate year may become less pronounced.

Figure 2 :

 2 Figure 2: Simulated time series. Attractor of sockeye salmon dynamics generated by computation of the timecontinuous model, showing oscillations with period 4 and the same pattern in the sequence of salmon abundance as in the empirical data. Only the biomasses at the end of each growth season are shown. a: deterministic model, b: dynamics with up to 50 percent fluctuations in the survivability of the sockeye salmon in the ocean (proportional to γ), see Discussion. Parameter values for these simulations are as explained in section 2.

Figure 3 :

 3 Figure 3: Bifurcation diagram of the biomass of the sockeye fry at the end of the season. The bifurcation parameter is the constant fraction of the zooplankton carrying capacity K 0 . For each value of K 0 , 200 consecutive data points are plotted.

Figure 4 :

 4 Figure 4: Eigenvalues of the matrix from Eq. (4), as obtained from a linear stability analysis of the time-continuous model. K 0 increases from 6.4 to 7.5 as the eigenvalues move outwards. The zoom shows that the eigenvalues close to the imaginary axis (+) are the first ones to cross the unit cycle.

Figure A1 :

 A1 Figure A1: Escapement data of for sockeye spawing stocks in the Fraser River basin. a: Quesnel Lake, b: Stuart Lake, c: Bowron Lake, d: Seton Lake, e: Chilko Lake.

Acknowledgements

This work was supported by the German Research Foundation (DFG) under contract number Dr300/7 and Br2315/9-1. Kim Hyatt, Jeremy Hume, and Carl Walters commented on an earlier version of this manuscript. Empirical data were kindly provided by Tracy Cone (escapement data) and by Kim Hyatt (age composition of the sockeye stocks). This collaboration started at the Complex Systems Workshop in Fairbanks, Alaska, 2007. The 2008 Cultus Lake Modelling Workshop, which had an important impact on this project, was hosted by Fisheries and Oceans Canada.