Modelling interaction patterns in a predator-prey system of two freshwater organisms in discrete time: an identified structural VAR approach

In ecology, the concept of predation describes interdependent patterns of having one species (called the predator) killing and consuming another (the prey). Specifying the so-called functional response of prey populations to predation is an important matter of debate which is typically addressed by means of continuous time models. Empirical regression or autoregression models applied to discrete predator-prey population data promise feasible steady state approximations of often complicated dynamic patterns of population growth and interaction. Ewing et al. (Ecol Econ 60:605–612, 2007) argue in favour of the informational content of so-called vector autoregressive models for the dynamic analysis of predator-prey systems. In this work we reconsider their analysis of dynamic interaction of two freshwater organisms, and design a structural model that allows to approximate the functional response in causal form. Results from an unrestricted structural model are in line with core axiomatic assumptions of predator-prey models. Conditional on population growth lagged up to three periods (i.e., 36 h), the semi-daily population growth of the prey Paramecium aurelia diminishes, on average, by 1.2 percentage points in response to an increase of the population growth of the predator Didinium nasutum by one percentage point.


Introduction
In the ecological and biological literature so-called predator-prey (PP) models have become an established framework to describe patterns of predation, i.e., the killing and consumption of one species (the prey) by another (the predator). The instantaneous rate of prey consumption per predator-the so-called functional response-is both an essential feature of PP interaction and an important dimension for structuring a broad variety of continuous time PP models (see Jost and Ellner 2000, for an overview). 1 The debate on the functional response (see Jost and Ellner 2000, for a comprehensive discussion) largely fluctuates around the functional specification and subsequent estimation of the so-called predation function in continuous time. A few studies have provided empirical assessments of PP models in discrete time. 2 Unlike continuous patterns of PP interaction, discrete time models can capture multiperiod dynamics in form of regressive or autoregressive patterns without explicit reference to an underlying theoretical model. Adding to this merit, it is worth to note that empirical sampling might take place at low frequencies such that important intergenerational dynamics could be subject to (implicit or unspecified) aggregation. In this context, Ewing et al. (2007) (henceforth, ERE07) are among the first to argue convincingly in favour of the informational content of vector autoregressive (VAR) models for the dynamic analysis in PP systems. 3 To materialize their claim, ERE07 provide an in-depth analysis of the semi-daily population data from the classic ciliate experiments of Veilleux (1979) in their digitalized form of Jost and Ellner (2000). In this work, we reconsider the analysis of ERE07. Apart from replication exercises, we complement their analysis with the quantification of instantaneous responses among prey and predator population growth rates in either direction.
ERE07 can be considered as the first scholars who illustrate dynamic model implications by means of impulse response analysis (or so-called innovation accounting). In particular, they adopt the generalized impulse response functions (generalized IRFs, GIRFs) of Pesaran and Shin (1998). Since GIRFs lack a strictly structural (or causal) interpretation (see, e.g. Kim 2013), however, ERE07 remain ultimately silent about a steady state assessment of the functional response. For example, GIRFs displayed in Fig. 2 of ERE07 suffer from the counter intuitive interpretation that positive surprises to the prey population invoke a significant deterioration of predator population growth. The structural model approach that we undertake in this work largely benefits from recent contributions to structural VAR (SVAR) analysis (see Kilian and Lütkepohl 2017, for an up to date textbook treatment of identification in SVARs). Specifically, we trace the stochastic variations in the PP system back to unique (i.e., identified) non-Gaussian independent components (Lancaster 1954;Comon 1994). The model implied impact effects of the independent components on population growth rates accord with the axiomatic assumptions of the so-called Kolmogorov model as the theoretical and intuitive foundation of many representatives of PP models (Freedman 1980). Conditional on the process history, by implication, an increase of the predator population by one percentage point reduces prey population growth by about 1.2 percentage points, on average.
The remainder of this work is organized as follows: The next Section sketches the data, highlights the (structural) identification problem and puts the independent component analysis (ICA) adopted in this work into the perspective of alternative approaches to identification in SVARs. The main empirical results are provided and discussed in Sect. 3. Section 4 concludes. The Appendix provides an explicit representation of the dependence coefficient of Bakirov et al. (2006) which we employ for the detection of independent components.

A structural model for predator-prey interaction
This Section develops the structural model of interest. On the one hand, we highlight stylized characteristics of the population growth rates under scrutiny and point to the structural VAR as a means to translate correlation estimates into causal linkages. On the other hand, we review briefly alternative approaches (theory-vs. data-based) to model identification, and argue in favour of ICA for the subsequent analysis.

The data and conditional correlations
The population data analysed in ERE07 have been determined in mixed ciliate experiments of B. Veilleux with two freshwater organisms, the prey Paramecium aurelia and the predator Didinium nasutum (Veilleux 1979). Specifically, with initial densities of 45 (Paramecium) and 15 (Didinium) individuals per milliliter (ml) ciliates were cultured in Petri dishes containing six ml of culture medium maintained at 27 C. To achieve sustained patterns of PP interaction, B. Veilleux elaborated on diverse scenarios combining controlled amounts of a bacterial nutrient for feeding the prey species (Cerophyl) with additions of Methyl Cellulose which governed thickness of the medium and thereby the movements of the ciliates. Nondestructive population counts of individuals per ml were made at the semi-daily frequency that accounts for intergenerational dynamics, since the number of fissions per day is larger than two for both species. While the fission rates might depend on the actual laboratory experiment (see Figures 2 and 3 of Veilleux 1979) further complexity of intergenerational dynamics arises from the fact that fission rates of the predator (Didinium) are by about 30% larger than those of the prey (Paramecium). To achieve sufficient accuracy of the population counts each quote is the average of eight independent countings. For intermediate choices of nutrition available to the prey species (i.e., cerophyl concentration) the experiments resulted in stable population patterns that are likely to reflect 'inherent properties of the system and not an artifact due to the experimental method' (Veilleux 1979, p. 796). 4 As in ERE07 the population growth rates of interest are defined as where N t and P t are individuals per ml of Paramecium (the prey) and Didinium (the predator), respectively. The observation index t indicates measurement intervals of 12 h. The total number of observations (N t ; P t ) is 52. Hence, the number of available growth rates (n t ; p t ) is 51. ERE07 provide ADF statistics in favour of stationarity of both growth rates. The graphical display of the growth rates in Fig. 1 shows two markedly outlying observations which might reflect the initialisation of the experimental design. Therefore, this work complements the analysis of full sample data in ERE07 (51 observations) with a robustness analysis for a subsample that excludes the first two recorded growth rates (49 observations). Table 1 documents mean growth rates and standard deviation (sd) statistics for the two samples. Moreover, Table 2 documents some conditional correlations obtained from ad-hoc regression models. Stylized bivariate regressions as well as richer dynamic structures obtain a common effect direction for both alternative choices of the dependent variable (n t or p t ). Nevertheless, the documented regression outcomes allow for an interesting insight into conditional correlations among both growth rates n t and p t . While the stylized bivariate regression with full sample information is characterized by positive coefficient estimates, these conditional correlations turn negative after conditioning on the set fn tÀi ; p tÀi ; i ¼ 1; 2; 3g as additional explanatory variables. Apparently, the model estimation by means of stylized bivariate regressions suffers from omitted dynamics. Put differently, dynamic effects appear as an important building block for the understanding of the PP system at hand. Corroborating established features of linear projections, however, the (dynamic) regression models are apparently not suitable to discover causal relationships by means of conditioning n t on p t (or p t on n t ). It is likely that with semi-daily sampling the populations of the species (and, hence, population growth rates) are subject to joint endogeneity. As a result, using contemporaneous population or growth rate data as explanatory variables will induce correlation with regression residuals and, hence, estimation bias. Bias free estimation approaches would deserve instrumental information or the a-priori knowledge of the causal structure. Both preconditions hardly apply in the present context of the analysis of PP interaction among freshwater organisms. 4 See Table 1 of Veilleux (1979) for a summary of (un)stable dynamic behaviours within Methyl Cellulose cultures. Based on simulated data from a continuous time model of PP interaction Jost and Ellner (2000) provide further support to consider the data as drawn from stable states of a PP model. The population counts are available from the web, http://cbi-toulouse.fr/images/upload/prslbdatasets.txt. As ERE07 we analyse the last series provided ('cerophyl concentration, CC = 0.5' ; digitalized from Figure 14c of Veilleux (1979)).
Seeing that cross equation linkages hide causal structures in ad-hoc (dynamic) regressions, the subsequent steps do not only take account of dynamic profiles but also aim at a regression design composed of unrelated equations. For this purpose we first develop the VAR model of ERE07 further towards a structural representation tracing model residuals back to orthogonal (and identified) shocks. Subsequently, the isolation of these shocks provides a view at PP linkages in a form of unrelated equations.

The VAR model and its structural form
Let the vector y t collect the growth rates defined in (1) as y t ¼ ðn t ; p t Þ 0 such that quotes on prey (predator) population growth are ordered first (second) in y t . To analyse competitive struggle in a PP model, the bivariate VAR model of order q (VAR(q)) reads as (see Lütkepohl 2007, for a textbook treatment of VARs). Taking account of unconditional effects m is a vector of intercepts, and e t is an uninformative model residual with mean zero and covariance R, i.e., e t $ ð0; RÞ. The model parameters (A i ; i ¼ 1; . . .; p; and R) and residuals (e t ) in (2) can be consistently estimated by means of OLS estimation. Conditioning on predetermined variables, however, the  (2) does not provide an explicit representation of the contemporaneous linkages that operate among the two species.
To address structural or contemporaneous relations, it has become a widespread strategy to regard the residuals in e t as the outcome of a linear and non-singular weighting scheme that applies to latent and orthogonal shocks collected in a vector n t . Rewriting the model residuals accordingly as e t ¼ Dn t obtains the SVAR counterpart of (2), Without loss of generality, the orthogonal shocks in n t are often assumed to have (marginal) variances of unity, i.e., Cov ½n t ¼ I with I denoting the identity matrix. The residual covariance R from (2) and the weighting parameters in D (see (3)) are related such that Highlighting the core identification problem, it is easy to see that the space of potential covariance factors D ¼ fDjDD 0 ¼ Rg is infinite. For instance, a parametrized space of covariance decompositions obtains from where G is a lower triangular Cholesky factor (R ¼ GG 0 ) and is a rotation matrix with rotation angle h. Accordingly, the covariance decomposition space could be defined as D ¼ fD h jD h ¼ GR h ; 0 h pg, and implied orthogonalized residuals are n ðhÞ

Identification
The identification problem in structural VAR analysis consists of the selection of one specific D matrix from D. Since all members of D align with the model outlines in (2) and (3), it is immediate that the identification problem in structural VARs cannot be solved without further external assumptions. In the case of a conditional Gaussian model, for instance, all rotations of n t (or equivalently of D) are informationally equivalent. Identifying external information might stem from two sources, (i) theoretical considerations about the variables in y t (i.e., about PP responses to the elements in n t ), or (ii) from stylized features of the distributions of e t or n t . We next encounter four alternative identification schemes in a stylized manner. For a detailed and up-to-date textbook treatment of identification in SVARs we refer the reader to Kilian and Lütkepohl (2017). Throughout we consider the case of bivariate SVARs. Extending the arguments to higher model dimensions is straightforward. Subsequently, we briefly discuss particular merits of theory-vs. data-based identification, and argue for ICA as a means to identify the PP SVAR model. Finally, we provide a short outline of the extraction of independent components.

Identification schemes
1. Recursive systems: The estimate of the covariance R comprises three single (co)variances. To identify the four structural parameters in D, Sims (1980) has suggested to impose a (lower) triangular structure on this matrix. With such a restriction, the unrestricted parameters in D can be quantified straightforwardly. However, this identification scheme comes with the imposition of a hierarchical structure of the model. By implication, the structural properties depend on the ordering of the variables in y t . 2. Identification by means of sign restrictions: As formalized in (3) typical elements of D quantify the impact effects of a unit shock n it ¼ 1; i ¼ 1; 2; on the observables in the system. While this effect is unknown in its exact metric form, it is often the case that theoretical considerations imply specific effect directions. In the econometric literature Faust (1998), Canova and De Nicolo (2002) and Uhlig (2005) have pioneered and advocated the identification of SVARs by means of such consensual effect directions. In the framework of PP models a theory guided shock labelling could benefit from theoretical knowledge of how the organisms under scrutiny react to surprising changes of living conditions (i.e., shocks). For instance, in real life scenarios such shocks might materialize in the form of the invasion of new species or changing environmental conditions. In this context it is worth observing that many modern representatives of PP models are based upon the Kolmogorov model -a minimal set of two differential equations coupled with axiomatic assumptions (Freedman 1980). Population changes as implied by the Kolmogorov model are and the axiomatic assumptions concern the so-called predation function g ( og oP \0) and, moreover, oh oN [ 0. Approximating the continuous time model in discrete time with dt ¼ 1 obtains n tþ1 % gðN t ; P t Þ and p tþ1 % hðN t ; P t Þ. Hence, a positive shock which primarily impacts on preys implies n tþ1 [ 0 (direct effect) and p tþ1 [ 0 (indirect effect channelled through h). Moreover, shocks which are directly beneficial to growth rates of the predator population are likely to impact adversely on the living conditions of the preys (indirect effect channelled through g). 5 Casting these considerations in the framework of the discrete time SVAR model obtains the following effect directions for typical Hence, a theory supported sign pattern for the structural matrix D is Providing a clear and unique effect pattern, the sign restrictions in D theo qualify for theory based model identification. Practically, a set of models identified by means of sign restrictions obtains from random sampling of rotation angles h; 0 h p, and storing all implied covariance decompositions that align with the sign pattern provided in (4). Once an analyst has sampled a sufficient number of admitted structural matrices, the set of identified models e D might be described by common statistical measures (e.g., mean, standard deviation, quantiles). 3. Identification through heteroskedasticity: The identification through heteroskedasticity as pioneered by Rigobon (2003) and developed further by, amongst others, Lanne and Lütkepohl (2008) and Bacchiocchi and Fanelli (2015) exploits information on observable heteroskedasticity in the residuals e t . Assume that an analyst has access to two distinct covariance estimates, denoted R ð1Þ and R ð2Þ , that she obtains from disjoint subsamples of the data. In the bivariate case this yields six single (co)variance estimates to solve the identification problem. The two covariances allow for a reparameterization as R ð1Þ ¼ DD 0 and R ð2Þ ¼ DWD 0 , where W is a diagonal matrix. Hence, the six estimated (co)variances can be mapped in a one-to-one manner to the six unknown parameters in D and W, thereby identifying D. It is evident, however, that unique identification by means of heteroskedasticity requires that the parameters in W must be different. In other words, the switch from the first to the second covariance regime must not affect the variances of the shocks proportionally. 4. Identification by means of independent components: Early results on the properties of linear combinations of (non)Gaussian random variables have motivated to obtain structural shocks as independent components. While Lancaster (1954) has stated that rotation invariance is unique to independent Gaussian random variables, results of Comon (1994) imply that the weighting scheme e t ¼ Dn t allows a unique recovery of the independent shocks in n t from e t if and only if at most one element of n t exhibits a Gaussian distribution. Taking advantage of these results, amongst others, Moneta et al. (2013), Gouriéroux et al. (2017) and Lanne et al. (2017) have developed identification schemes for SVARs which differ with respect to the space of admitted D matrices or the distributional assumptions made for the structural shocks in n t . Non-parametric approaches to ICA have been suggested or applied, for instance, by Matteson and Tsay (2017) and Herwartz (2019). Following their lines of reasoning an estimator of D could be obtained from solving the minimization problem by means of a general non-parametric test of the null hypothesis of joint independence.

A comparative discussion of identification schemes
To support the choice for a particular identification scheme in empirical practice it is helpful to first contrast theory-vs. data based identification. While theory-based identification ensures that obtained structural shocks have meaningful impact properties, such schemes are also prone to intermingling model assumptions and conclusions (Uhlig 2005). In contrast, data-based identification typically provides sufficient external information such that otherwise overidentifying restrictions could be subjected to testing. At the downside, however, structural shocks obtained from these identification schemes do not necessarily allow for a sound theoretical interpretation. In this work we opt for a data-based approach to identification, and investigate subsequently if the retrieved shocks allow for a theory-conform interpretation in the spirit of the theoretical sign pattern of D theo in (4). For this purpose it is important to notice that an identified matrix D is unique only up to column signs and the column ordering. With d i denoting the i-th column of D it Hence, the matrices ½d 1 ; d 2 , ½Àd 1 ; d 2 or ½Àd 2 ; d 1 are all in line with the covariance restriction R ¼ DD 0 .
The application of data-based identification schemes deserves either that the VAR residuals are subject to informative covariance changes or non-Gaussian. While we obtain some evidence in favour of downward shifting variances of VAR residuals (see also Fig. 1), the null hypothesis of proportional variance changes cannot be rejected with 5% significance. In contrast, the evidence against normally distributed VAR residuals is very strong such that we choose ICA for identification. The following paragraph is explicit on the independence statistic used to select D from D.

Independent component analysis
Solving the minimization problem in (5) requires a statistical quantification of joint dependence. To detect unique independent components we follow the approach of Herwartz (2019) and employ the dependence coefficient introduced by Bakirov et al. (2006). The dependence coefficient, denoted C, provides a ffiffiffi ffi T p consistent estimator of the population distance of the joint characteristic function of two random variables (in our case n 1t and n 2t ) and the product of both marginal characteristic functions (see the Appendix for an explicit exposition). Similar to absolute estimates of linear correlations, C is bounded between zero and unity which indicate the limiting scenarios of independence and complete dependence, respectively. For testing the null hypotheses of independence, TC 2 is a bounded random variable. Since critical values of TC 2 are difficult to obtain analytically, Bakirov et al. (2006) suggest a permutation bootstrap approach for the assessment of significance of the test statistic. 6 To determine independent components practically, one might either minimize C or maximize the p-value of TC 2 . In light of simulation results of Herwartz (2019) we opt for the latter and determine the structural matrix where C h is determined from the sample fn ðhÞ t ¼ D À1 h e t g T t¼1 . 7 Practically, the estimatesĥ obtain from a grid search over alternative rotation angles h ¼ jp 2 =50; j ¼ 1; 2; 3; . . .; 50. Extending the grid to p 2 \h p results in an informationally equivalent second solution for D h with exchanged column positions. 6 For all computational purposes that involve C we use the R package 'energy' (Rizzo and Székely 2017). Bootstrap p-values for TC 2 are based on R ¼ 249 permutations. It turns out that subjecting a given sample to independence testing obtains some variation in the implied p-values for a given statistic TC 2 . To assess the significance (and p-values) of TC 2 robustly, all respective resampling exercises build upon average p-values from 49 alternative initializations of the permutation bootstrap. Results obtained from using just one initialization for the permutation bootstrap are qualitatively identical, however. 7 Solving an estimation problem by maximizing the p-value of a suitable test statistic (i.e., minimizing the evidence against the null hypothesis) has become prominent as so-called Hodges-Lehmann estimation (Hodges and Lehmann 1963).

Marginal effects
The model in (3) is structural in the sense that observable regression residuals are traced back to orthogonal shocks. As such the model is not explicit on the implied causal relations between the jointly endogenous variables in y t ¼ ðn t ; p t Þ 0 . Multiplying (3) from the left with D À1 obtains a structural model with unrelated residuals and a parametric representation of the conditional linkages among the elements in y t , where . . .; q. Let X tÀ1 denote an information set that comprises the process information up to time t À 1. Moreover, let d ðijÞ denote a typical element of D À1 . Conditional on X tÀ1 , the elements in D À1 describe the marginal causal effects. Since we obtain after normalization from (6) where x ðnÞ tÀ1 and x ðpÞ tÀ1 are observable conditional on the VAR parameters and X tÀ1 . It is worth noticing that by virtue of their definition in (7) the marginal effects Àd ð21Þ =d ð22Þ obey a similar interpretation in discrete time as the predation function g(N, P) in the continuous time Kolmogorov model.

Impulse response functions
Owing to the rich and complex dynamic parametrization of VAR models, IRFs have become a common practice to illustrate model implied dynamics. To formally derive IRFs as in Chapter 2 of Lütkepohl (2007) the VAR in (2) can be represented in so-called moving average form as where L is the backshift operator such that Ly t ¼ y tÀ1 and l ¼ A À1 ðLÞm. The moving average parameter matrices W i ; i ¼ 0; 1; 2; 3; . . ., can be obtained recursively from nonlinear transformations of the autoregressive parameter matrices A i ; i ¼ 1; . . .; q. Specifically, The moving average representation is of core importance for purposes of impulse response analysis. For instance, the typical elements of the matrices characterize the effect of isolated unit shocks in n jt on y i;tþh , the i-th time series variables at horizon h. Obviously, impulse response schemes are specific to the selection of the structural matrix D.
To avoid strong identifying assumptions on the covariance decomposition, the innovation accounting analysis in ERE07 does not rely on impulse responses to orthogonalized shocks. Rather ERE07 make use of the so-called GIRFs of Pesaran and Shin (1998). Formally, GIRFs read as where e i is a zero vector with a unit element in its i-th position, and r ii is the square root of the i-th diagonal element of the residual covariance R. The qualification of the profiles in (10) as 'generalized' IRFs comes from the fact that these statistics are invariant to the ordering of the variables in the system. As argued in Kim (2013), however, GIRFs could be interpreted as an aggregation of impulse responses from VARs with distinct variable orderings. By implication, GIRFs lack effect consistency under the typical case of a non-diagonal residual covariance R.

Structural analysis of predator-prey interaction
In this Section we first highlight that the VAR model of ERE07 is characterized by significant deviations from conditional normality and, hence, that independence scores might be used for characterizing ad-hoc identification schemes in terms of independence of orthogonalized shocks. Subsequently, we pursue the detection of independent components and provide an interpretation of the identified shocks. Furthermore, the Section includes a comparative description of the GIRFs of ERE07 and identified IRFs. Finally, we employ the identified model to assess the contemporaneous causal linkages among the growth patterns of the two species.

Non-Gaussianity and recursive model structures
ERE07 argue in favour of a VAR model with q ¼ 3, as it minimizes the BIC and obtains serially uncorrelated estimates of e t . After fitting a VAR(3) model to the bivariate vector of population growth rates, the residual covariance matrices estimated from the full and reduced sample read, respectively, as While the removal of initial outlying observations reduces markedly the residual variances of both variables, the residual correlations are significantly negative, namely -0.477 and -0.458 for the full and the reduced sample, respectively. To uncover the causal relationships behind these empirical correlations by means of unique independent components, it is essential that the data are non-Gaussian. Table 3 provides outcomes from normality tests applied to VAR residuals. Both the Shapiro-Wilk (SW) and the Jarque-Bera (JB) test statistics for the null hypothesis of normally distributed marginal residuals provide strong evidence against the normality assumption. JB summary statistics indicate non-Gaussianity also for the joint distribution. Given significant deviations from the Gaussian model, it is tempting to order alternative ad-hoc orthogonalization schemes with respect to their potential to align with notions of independent (rather than merely orthogonal) shocks. Cholesky factors of the residual covariance depend on the ordering of the variables in y t (see the discussion in Section 2.3.1). Alternatively, for a given variable ordering one might distinguish a lower (denoted D l ) and an upper triangular Cholesky factor (D u ). The left hand side and middle panels of Table 4 document estimated alternative covariance decomposition factors. In terms of joint independence the orthogonalized shocks implied by alternative hierarchical model structures allow a clear distinction. Considering residuals e 1t (prey growth) to originate in isolation through orthogonalized shocks (as implied by D l ) is in line with the null hypothesis of independence at conventional significance levels. While one could argue that accepting a null hypothesis as such does not provide a strong inferential result, it is worth to highlight that the data have no possibility to object against an The table displays results from Shapiro-Wilk (SW) and Jarque-Bera (JB) tests. The latter also include diagnostics for the joint distribution and a separation for testing skewness and excess kurtosis separately.
p-values are multiplied with 100 orthogonalization of a model with prey residual growth ordered first. Unlike this hierarchical specification, assuming residuals e 2t (predator growth) to originate in isolation through orthogonalized shocks (as implied by D u ) obtains a rejection of the independence hypothesis with 5% (10%) significance for the full (reduced) sample. In sum, statistical identification is an informative means to establish diagnostic evidence in favour of the actual variable ordering y t ¼ ðn t ; p t Þ 0 over the alternative y t ¼ ðp t ; n t Þ 0 . 8

Identified independent components
The right hand side panel of Table 4 documents structural parameters estimates b D joint with bootstrap-based t-ratios. 9 For both samples the detection of shocks with minimum dependence results in solutions which are close to the orthogonalized shocks obtained from the Cholesky covariance factor. The maximization of model implied p-values obtains estimates of about 60%. While one should be careful in interpreting supremum p-values from multiple testing in the usual way, it seems that the shocks implied by b D show some weaker dependence in comparison with the shocks retrieved after imposing a hierarchical model structure (D l ). However, this improvement is marginal and might lack significance.
The estimator b D provides evidence for two well identified shocks, since the sign patterns of estimated effect directions are distinct. While the first shock invokes opposite impacts on the two population growth rates, the second shock exerts a positive impact effect on both the growth rate of the predator and the growth rate of the prey population. Supporting the restriction implied by the Cholesky factor D l , however, the upper right element in b D lacks positiveness with conventional significance.
The estimates in b D are the solution of a statistical optimization problem. As resulting from a purely data-driven approach, however, it is not clear if the implied structural shocks fn t ¼ b D À1 e t g T t¼1 conform with any theoretical interpretation. In data-based SVAR identification the issue of shock labelling has become an important step of the analysis. In the form provided in Table 4 both structural 8 See also the simulation results in Herwartz (2019) on the informational content of the independence coefficient C to rank causal hierarchies implied by alternative variable orderings. 9 For inferential assessments of estimates in b D (and (G)IRFs) we use a fixed-design wild bootstrap scheme in the spirit of Gonçalves and Kilian (2004). At each bootstrap replication pseudo samples are composed of y Ã t ¼m þÂ 1 y tÀ1 þÂ 2 y tÀ2 þ Á Á Á þÂ p y tÀp þ e Ã t ; e Ã t ¼ê t w t ; t ¼ 1; 2; . . .; T: In (12)m andÂ j ; j ¼ 1; . . .; q, are OLS parameter estimates retrieved from the data. Independently of the data, the scalar random variable w t exhibits a Rademacher distribution ( Prob ½w ¼ 1 ¼ Prob ½w ¼ À1 ¼ 0:5). The inferential analysis relies on 999 bootstrap replications. Brüggemann et al. (2016) have shown that wild bootstrap variants fail to mimic the joint distribution of estimators of slope (m; A i ; i ¼ 1; 2; . . .; q) and (co)variance parameters (R ¼ DD 0 ) asymptotically. In finite samples, however, a consistent moving block bootstrap and the wild bootstrap likely obtain qualitatively identical results (Brüggemann et al. 2016). matrices D l and b D point to the prominent role of the first shock (i.e., isolated shocks entering the model through prey growth rate residuals). As implied by the documented sign pattern, this shock improves prey growth and diminishes growth rates of the predator population. In a 'material' sense, such theoretical properties of 'positive' shocks might lack intuition. Therefore, we reorganize the columns and sign patterns of b D to establish features of the shocks that align with D theo in (4). Specifically, we reverse the column ordering and subsequently multiply the elements of the second column with minus unity. This obtains the structural matrix estimates in labelled form for the full sample and the reduced sample, respectively. When it comes to describe actual shocks behind the empirical growth rates analysed in this work, the controlled conditions of the experiments described in Veilleux (1979) leave little room to assign a theoretical interpretation to the labelled structural shocks. For instance, the experiments control the amount of food that was available to the prey species (i.e., the cerophyl concentration). However, potential 'shocks' which improve the living conditions of the prey species could occur in the forms of (i) favourable stochastic fluctuations of the nutrient (or its quality) around the controlled mean level, or (ii) an advantageous spatial distribution of the nutrient which minimizes search time/cost to the prey population. Improving growth rates of the Paramecium (prey) population, such shocks also improve the ingestion conditions for Didinium (predator). Therefore such shocks can be seen to invoke impact responses as described in the first columns of the structural matrices in (13) for the identified n 1t . For arguing in favour of potential positive shocks to the predator population it is worth noticing that the predator Didinium is unable to detect its prey Paramecium visually or in a chemotactic manner. Hence, it cannot actively move towards its prey, and it necessitates a random direct contact to induce the capture and ingestion mechanism. Therefore particular states of species distribution that are favourable for predator growth include the absence of predator clustering or the presence of prey clustering. Establishing favourable conditions for predator growth, such 'shocks' are most likely detrimental for prey growth, and point to the impact effects in the second columns of the structural matrices in (13) and the properties of the identified n 2t .

Impulse response analysis
We next describe the implications of the VAR model for the dynamic and contemporaneous interplay of the two species. For this purpose we discuss a replication of the GIRF estimates of ERE07 which is displayed in Fig. 2. As a complement to these results, we discuss effects implied by the identified SVAR in the form of IRFs to orthogonalised and independent shocks shown in Fig. 3.

Generalized impulse response functions reconsidered
While point estimates shown in Fig. 2 fully accord with the GIRFs of ERE07, the displayed confidence intervals with 95% coverage differ slightly in terms of shape. 10 The displayed functional forms indicate the stationarity of the system of  (10)) of Paramecium (prey, n tþh ) and Didinium (predator) population growth rates (p tþh ) to unit shocks. The time index t operates at the semidaily frequency. Point estimates are equivalent to those shown in Figure 2 of ERE07. Inferential assessments rely on a fixed-design wild bootstrap scheme in the spirit of Gonçalves and Kilian (2004) (see footnote 9 for details). Dashed curves indicate point-wise 95% bootstrap confidence bands (i.e., the 2.5% and 97.5% quantiles of the bootstrap distribution). Displayed GIRFs are for the full sample (FS). Results for the reduced sample (RS) are almost identical population growth rates as the effects of initial impulses die out over an horizon of approximately seven days. As highlighted in ERE07, features of both persistence and cyclicality characterize the dynamic system. Seemingly, prey growth rates are in the lead of predator growth rates.
Positive impulses in the growth rates of the predator or prey population invoke markedly negative impact effects in the growth rates of the other population. This outcome is a direct consequence of the negative contemporaneous correlation characterizing the residual covariances documented in (11). Hence, the GIRFs lack a strictly structural interpretation. In the present case this lack shows up in an important theoretical model deficiency. While it is most intuitive to see that positive shocks to the predator population growth are -on impact -detrimental for the growth rate of the prey population, the reverse transmission channel lacks such an intuition. In fact, the model implication that positive shocks to prey population growth exert a negative impact on the growth rate of the predator population Fig. 3 Impulse responses of Paramecium (prey) and Didinium (predator) population growth rates to identified orthogonal and independent unit shocks (see (9)). For further notes see Fig. 2 contradicts basic assumptions of a stylized Kolmogorov PP model. This result is supportive for the critical perspective on GIRFs provided by Kim (2013).

Identified impulse responses to structural shocks
Continuing the structural analysis we discuss the dynamic response patterns of prey and predator population growth rates to stylized (orthogonalised and independent) unit shocks. With regard to the first shock the upper left panel of Fig. 3 indicates that its positive impact on prey growth turns significantly negative within the first 24 hours. Then, cyclical effects imply a return to positive rates of population growth which lack 10% significance. The response of the prey population to this shock leads the response of the predator population growth. The latter achieves its minimum value after 48 hours and the maximum of recovery also follows the maximum prey recovery with a delay of 24 hours.
The second shock displayed in the right hand side panels of Fig. 3 impacts positively on predator growth and shows a markedly diminishing effect on prey growth. Due to the deterioration of the availability of nutrition, the initially positive response of predator growth turns negative within 24 hours to allow an almost significant recovery of the prey population followed again by positive dynamic responses of predator population growth.

Marginal effects
As argued in Sect. 3.4 the structural model and -in particular -its reformulation in (6) allows to recover the contemporaneous effect of growth rates in either species on the remaining one conditional on the history of the process. Structural model estimates of inverse (and normalized) b D À1 matrices are D À1 ðFSÞ ¼ for the full sample and the reduced sample respectively. The estimated structural model implies the contemporaneous relationships between prey and predator population as documented in Table 5. In the steady state an increase of the predator population growth by one percentage point invokes a reduction of the prey Bootstrap interquartile ranges of parameter estimates in parentheses. The causal patterns of interest are indicated with '! ' . Causal effect estimates p t ! n t (n t ! p t ) obtain from inverted structural estimates after normalization as Àd ð12Þ =d ð11Þ (Àd ð21Þ =d ð22Þ ), where d ðijÞ is the ij-element of D À1 population growth by about 1.2 percentage points. Bootstrapping the structural model obtains an interquartile range of this relation which is between 1.06 and 1.48 (1.02 and 1.39) percentage points for the full (reduced) sample. Regarding the reversed directional impact, an increase in prey growth by one percentage point increases, on average, predator growth by 5.1% (6.8%) percentage points for the full (reduced) sample. Although these estimates are representative for the steady state, it is worth noticing that in directional terms both contemporaneous causal effects are in line with the axiomatic assumptions of the Kolmogorov model.

Conclusions
In the ecological literature predator-prey (PP) models have become a successful approach to understand species' interactions. Ewing et al. (2007) have argued convincingly in favour of the informational content of vector autoregressive models for the dynamic analysis in PP systems conditional on samples drawn (with low frequency) at discrete time. In this paper we complement their analysis of a prominent laboratory time series of population counts of two freshwater organisms (Veilleux 1979). Generalised impulse response patterns (Pesaran and Shin 1998) displayed in the benchmark study of Ewing et al. (2007) are not fully in line with the axiomatic foundations of PP models as stated in the Kolmogorov model (Freedman 1980). This caveat can be traced back to the fact that generalised impulse response functions lack a strictly structural interpretation.
We take advantage of recent contributions to the literature on identification in structural VAR models to provide insights into both the dynamic and the contemporaneous interplay of the two freshwater organisms. In specific, our analysis builds upon the uniqueness of independent shocks which exhibit a non-Gaussian distribution (Lancaster 1954;Comon 1994). Independent components give rise to a structural model which is well identified and in line with the basic axiomatic foundations of PP models. The statistical approach to identification allows to test the otherwise just identifying assumption of a hierarchical model structure. If the hierarchical model develops from regarding surprises to prey growth as rescaled structural shocks, the model is in line with an assumption of independent shocks. Estimating the structural parameters unrestrictedly provides further support for the hierarchical structure. With regard to the functional response in PP models our steady state approximation implies that conditional on a history of three half days, the prey population growth 'instantaneously' shrinks by about 1.2 percentage points in response to an increase of the predator population by one percentage point.
Adopting statistical identification schemes to the laboratory data of Veilleux (1979) shows the potential of structural perspectives on PP models of competition, coordination and predation. As a matter of fact, the laboratory design of data generation limits the material interpretation of the identified shocks. An interesting avenue for future research is the adaptation of structural VARs to real life data in ecology and other disciplines (e.g., the prominent linkage among 'snowshoe hare' and 'Canada lynx ' Stenseth et al. 1997). In such a context interesting ecological issues could be addressed by means of SVARs that are identified in a statistical and/ or theory based manner. For instance, the effects of controlled releases of the predator species into the wild could be assessed by means of structural impulse responses. Similarly, to understand the possible effects of weakened protection rules for preys (and, hence, more hunting) structural impulse responses promise valuable information.
Appendix: The dependence coefficient Let n st ¼ ðn 1s ; n 2t Þ 0 ; t ¼ 1; 2; . . .; T denote a sample of bivariate tuples consisting of two random variables n 1t and n 2t . Then, the dependence coefficient of Bakirov et al. (2006) is where jjn st À n pq jj; jn 2s À n 2t j and w 5 ¼ 1 jjn pp À n st jj; and j Á j and jj Á jj are the Euclidean norms in R and R 2 , respectively. It is worth noticing that in the empirical application pursued here, estimated quantitiesn t enter C and not directly observed random variables. However, from the ffiffiffi ffi T p consistency of estimated VAR residualsê t it is easy to show that observation specific estimation errors ofn t vanish on average when determining the sample moments w i ; i ¼ 1; 2; . . .; 5. Unlike rank-based dependence diagnostics (as, e.g., Hoeffding's D) the dependence coefficient C is informationally efficient, since it relies on the original random variables in n t in their metric form.