Spatio-Temporal Analysis of Game Harvests in Sweden

The benefits and costs of wildlife are contingent on the spatial overlap of animal populations with economic and recreational human activities. By using a production function approach with dynamic spatial panel data models, we analyze the effects of human hunting and carnivore predation pressure on the value of ungulate game harvests. The results show evidence of dynamic spatial dependence in the harvests of roe deer and wild boar, but not in those of moose, which is likely explained by the presence of harvesting quotas for the latter. Results suggest the impact of lynx on roe deer harvesting values is reduced by 75% when spatial effects are taken into account. The spatial analysis confirms that policymakers’ aim to reduce wild boar populations through increased hunting has been successful, an effect that was only visible when considering spatial effects.


Introduction
Conflicts between wildlife and economic activities is a worldwide phenomenon. Some of the major reasons are large carnivores' predation of livestock and wild game, and ungulates' browsing damage on agricultural and forest crops, creating conflicts with pastoralists, hunters, and farmers, respectively (Gren et al. 2018(Gren et al. , 2020Treves and Karanth 2003;Shwiff et al. 2013). Conservation of wildlife is a major aim in most parts of the world, and involves benefits and costs to society that vary depending on the spatial distribution of the animal species populations, landscape, and the human activities that are affected (Boman et al. 2003;Treves et al. 2004;Van Eeden et al. 2018;Widman and Elofsson 2018). Different perceptions of the magnitude of costs and benefits of wildlife among different stakeholder groups can often fuel other urban-rural and political conflicts (Krange et al. 2022).
The role of spatial spillovers for the benefits and costs from policies for wildlife conservation and management are relatively little studied. Boman et al. (2003) use a static economic optimization model for wolf management, taking into account wolves' withinspecies density-dependent dispersal. Aadland et al. (2016) develop a spatial-dynamic predator-prey economic model on wolf and elk management, assuming that prey dispersal is a function of predation risk, while Aadland et al. (2015) model the optimal management of a pest, the mountain pine beetle, for which the dispersal is assumed to be driven by the density of their prey, in this case trees. Rollins and Briggs (1996) analyze the optimal design of contracts for wildlife damage prevention and compensation, when efforts to prevent wildlife damage at a given location risks leading to increased damage in neighbor locations. A couple of studies show that economic activities could relocate in response to wildlife density, e.g., farming is shown to relocate to regions less prone to wildlife damage (Hansen et al. 2019;Mmopelwa and Mpolokeng 2008), and hunters incur additional costs for travelling to locations with a higher abundance of game (Knoche and Lupi 2013;Sarker and Surry 1998). Also, Mensah and Elofsson (2017) identify the presence of spatial spillovers in hunting lease prices. With applications to fishery, a couple of studies elaborate theoretically on the economically optimal management of spatial-dynamic systems (see, e.g., Wilen 2005, 2007;Smith et al. 2009), accounting for within-species density-dependent dispersal as well as the spatial reallocation of the profit maximizing harvesting effort. Also, Abernethy et al. (2007) study empirically whether the spatial pattern of artisanal fishery is consistent with the predictions of the Ideal Free Distribution (IFD), requiring that the average profit per hour of fishing should be equal among similar fishers. Results showed this was not the case, and neither did the fishing pressure increase with resource availability.
Common to the mentioned studies is that species' dispersal is defined by an exogenously given parameter and a single driver. This contrasts with the ecological literature, which shows that the abundance and dispersal of a species is determined by multiple drivers including, e.g., the availability of food, the landscape structure (Saïd and Servanty 2005;Maier et al. 2005;Acevedo et al. 2014), and behavioral responses related to predator-prey dynamics (Houston et al. 1993;Kjellander et al. 2004;Lima and Zollner 1996;Hebblewhite and Merrill 2007;Bergerud et al. 1984) and human-wildlife interactions (Cromsigt et al. 2013;Thurfjell et al. 2017). Also, none of the studies acknowledges the potential impacts of human activities, such as hunting and wildlife policies (e.g., local harvest quotas and zoning policies for conservation of large carnivores), on wild species' dispersal propensity and the resulting spatial distribution of harvests. Thus, economic studies have neither explored the relative importance of different spatial-dynamic processes for the economic outcome, nor taken into consideration the multiple natural and economic drivers of these processes. This is a problem because ignoring these factors could lead to biased results. Without knowledge on these processes and their economic impacts, policymakers will not be able to foresee the economic and ecological consequences of policy decisions related to conservation of large carnivores and harvesting of wildlife, and the distribution of the effects across space.
The purpose of this paper is to analyze the effects of human hunting and carnivore predation pressure 1 on game harvests, acknowledging the spatial and temporal dimensions that arise due to ungulates' avoidance of predation and hunting risk, and carnivores' and hunters' seeking game. To analyze this issue, we apply a production function approach, first developing a dynamic bioeconomic model accounting for the linkages between harvests, human hunting effort, and carnivore abundance. We subsequently implement dynamic spatial econometric methods to analyze how carnivore and hunter abundance influence harvests of the three most harvested game species in Sweden: Roe deer (Capreolus capreolus), wild boar (Sus scrofa), and moose (Alces alces). The three species are subject to different policies: for moose, there are local harvesting quotas, roe deer harvests levels are unregulated, and wild boar hunting is encouraged by permitting more efficient methods and gear. Moreover, three large carnivores play an important role for the studied game populations: Wolves (Canis lupus) predate on both moose and roe deer, the omnivorous brown bear (Ursus arctos) predates on moose calves, and the relatively smaller lynx (Lynx lynx) predates on roe deer. We expect to find difference across species that reflect the differences in harvesting policy and in predator-prey relationships. The methodological approach in our study is new, as there is no economic study, or unbeknownst to us, that incorporates both dynamic and spatial econometric techniques when assessing wildlife management and conservation based on empirical data, aiming to quantify spatial spillover effects within an ecological-economic framework, and the associated impact on benefits and costs associated to wildlife management. Our findings indicate that one additional lynx individual can decrease the harvest of roe deer in the short run by 2.8 within the same county and by 2 in neighboring counties. Similarly, one added hunter reduces boar harvest by 0.35 within a county, while an additional wolf and an additional bear decrease the harvest of moose by about 15 and 4 respectively. Such marginal impacts are informative to policymakers that aim to achieve a balance between benefits and costs in carnivore conservation policies.
The remainder of the paper is organized as follows: Sect. 2 presents an empirical context of the hunting in Sweden, the game species and their large predators. Section 3 lays out a bioeconomic model that develops the theoretical benchmark of our analysis. Section 4 presents the econometric methodology used to investigate the spatial and time dynamics of the studied system, the econometric model selection criteria, and the method to calculate the marginal effects for the significant independent variables. Section 5 explains the data, and Sect. 6 presents the results. Section 7 provides a discussion and concludes.

Hunting in Sweden
Recreation, game meat supply, and wildlife management are the main motivations for hunting in Sweden. There are about 300,000 hunters in the country, whose hunting activities are directly or indirectly associated with game population control and ecosystem management practices (Lindqvist et al. 2014). Most hunters hunt in the neighborhood of their home. Landowners own the right to hunt on their land, and hunters can access such hunting land by paying a lease to the landowner. This payment of the lease price transfers the hunting rights to the hunters and allows them to keep the game harvested within the property. Typically, the lease entitles the hunter to hunt all species defined as eligible according to the Hunting Ordinance on the land in question. Lease contracts are usually on a yearly basis and many hunters lease the same ground over long time periods, while others hunt when they are invited as guests or purchase short term leases. The roe deer, the wild boar, and the moose are three of the most popular hunted game animals in Sweden. About 104,000 roe deer, 98,000 wild boar, and 83,000 moose were harvested in 2015 (www. viltd ata. se).

Ungulate Game Species
The roe deer is a medium-sized ungulate spread nationwide with lower population densities in the northernmost areas of the country. Harvests levels are unregulated, and in most regions the hunting season occur during autumn and winter, spanning from 4 to 5.5 months depending on sex and age of the animal (Elofsson and Häggmark 2021). Roe deer can inflict browsing damage on valuable trees and agricultural crops, but the magnitude of the damage is small compared to that for moose and wild boar (Bergquist and Örlander 1998;Elofsson and Häggmark 2021).
The wild boar is an ungulate game with high dispersal capacity but low survival rates of the offspring under severe winter conditions. Boar populations are absent in northern Sweden (Engelmann et al. 2016), but widely present in the middle and the south of the country. The wild boar can cause considerable burrowing damage to agricultural crops (Gren et al. 2020;Mensah and Elofsson 2017), and hunting of boar yearlings is permitted throughout the year, while adults can be hunted from April to February. In spite of this, population numbers are increasing and since 2018, more efficient hunting equipment such as variable lighting, light amplifiers and reticles with thermal imagery is permitted with an aim to control spread and abundance of boar.
The moose is a large ungulate widely present in most parts of Sweden. High moose populations provide large hunting benefits, 2 but increase the number of moose-vehicle collisions, and produce negative effects to forest crops due to browsing damages on valuable tree species (Angelstam et al. 2000;Seiler 2004). As a result, hunting moose is regulated through decisions at the regional (county) level, and yearly harvest rates are defined by the County Administrative Boards (CABs) for so-called Moose Management Areas (MMAs), which are territorial divisions confined by natural barriers and established to balance the interest of hunters with other stakeholders (Mensah et al. 2019, Franklin et al. 2020). There are several MMAs in each county, typically ranging between 50,000 and 100,000 hectares (Sandström et al. 2013). The hunting season is restricted to the autumn and winter.

Large Carnivore Species
Sweden is home of three large terrestrial carnivores: the wolf, the brown bear and the lynx. Wolves are mostly present in central Sweden and mainly preys on moose, but also on roe deer (Sand et al. 2016). In other countries, occasional predation on wild boar is observed 1 3 (Mori et al. 2017), but this is extremely rare in Sweden, potentially because dispersal areas hardly overlap. Brown bears are mostly found in central and northern Sweden. They are omnivorous, but frequently prey on moose and reindeer 3 calves in spring after hibernation (Dahle et al. 1998). Lynx, which are smaller than the two foregoing, are present all over the country, in particular in central Sweden where roe deer are abundant. The lynx mainly preys on roe deer, but can also take other small game such as hares (Andrén and Liberg 2015).
The three large carnivores are all protected under the European Union's Habitats Directive 91/43/EEC of February 2007, and hunting is prohibited unless the population substantially exceeds favorable conservation status. 4 Within our study period, the Swedish Environmental Protection Agency (SEPA) and the CABs permitted license hunting of lynx and brown bear, but not wolf. On the other hand, wolves are subject to strict population control in the northern part of the country to avoid conflicts with reindeer herding.

Theoretical Model
In the following, we develop a theoretical model of time dynamics in game harvesting. This model serves as a basis for our approach to time dynamics in the econometric analysis, which is described in Sect. 4. The development of the game population is determined by the size of the game population P , hunting effort E , and habitat conditions L , where the latter expresses habitat or landscape features that affect the survival and reproduction of the species such as predation risk or feed availability. Based on Barbier (2007), the change in game stock, P t − P t−1 , is assumed to be given by the biological growth in time period t − 1 , F(P t−1 , L t−1 ) , net of harvesting in time period t − 1 , h(P t−1 , E t−1 ): Following Foley (2010), we assume a biological growth function with facultative habitat use, i.e., where the habitat conditions L affect the game species population growth, but the presence or absence of such conditions do not lead to extinction of the game population (Auster 2005): where r is the intrinsic growth rate, and K is the carrying capacity, which relates to the maximum population without hunting and in the absence of the habitat conditions L (Clark 1990). In our empirical application below, we assume that the habitat conditions L are associated with carnivore abundance. Accordingly, is a (negative) parameter expressing the sensitivity of ungulate game to the presence of different carnivores, due to the impact of the carnivores on prey mortality and on prey's habitat selection. 5 Following Boman (1995), Bostedt and Grahn (2008), Zabel et al. (2014) and Skonhoft (2006), we assume that depredation has no feedback effects on the size of the carnivore populations. This assumption is motivated given that carnivore population sizes in Sweden are determined by political decisions.
Further, we adopt a Schaefer harvest-effort relation (Schaefer 1954): where q is the catchability constant, representing the hunting technology. Thus, harvest rates are increasing in catchability, game population size, and hunting effort in the same time period. Hunting effort is assumed to be exogenous, which is motivated by the focus on spatial spillovers and model tractability, in combination with the lack of suitable data for empirical modelling of hunters' complex choice of effort. 6 Substituting (2) and (3) into (1), we obtain: Dividing by P t−1 we get: Since there are no data on game populations P, we cannot estimate Eq. (5) directly. Instead, we follow Schnute (1977) by defining the harvest (or catch) per unit effort, c t , as: Plugging this expression into (5), we obtain: Rewriting Eq. (6), we obtain the following equation: which can be redefined and estimated using the following regression model: (8) provides a workable theoretical model, which explains that the rate of change of the harvest per unit effort ( S t ) is a linear and dynamic function of the four terms on the right hand side. First, the intercept 0 is expected to be positive given that the carrying capacity K and the intrinsic growth rate r are both positive numbers. Second, the term 1 L t−1 describes the impact of carnivores on the game harvest, hence the 1 coefficient is inferred to be negative because carnivores reduce the game population. Third, the term 2 c t−1 designates the effect of the lagged harvest per unit effort. Because q and r are both positive, the theoretical model predicts 2 to be negative. Fourth, the term 3 E t−1 describes the impact of hunting effort on the harvested game and given that q is positive, 3 is inferred to be negative.
The theoretical model expressed in Eq. (7) or (8) could also be interpreted using a level form mathematical formulation. 7 By solving for c t in Eq. (7), it can be shown that the harvest per unit effort equates to the following level form model specification: Equation (8*) can be interpreted in a very similar fashion as Eq. (7) because it is influenced by the same explanatory variables. However, the dependent variable is expressed in levels ( c t ) as opposed to Eq. (7) which is expressed in change rates . This implies that Eq. (8*) is a function of interaction and quadratic terms, hence non-linear. This is an important justification to adopt our empirical strategy, and it is carefully addressed in Sect. 4 hereafter.

Methodology and Empirical Strategy
Although Eq. (7) or its level form equivalent (8*) represents an interesting theoretical model explaining the per unit harvest effort, its econometric implementation raises several issues and questions which need to be resolved. These issues are now first discussed before we proceed with the empirical analysis. A preeminent problem stems from the fact that Eq. (7) does not account for spatial interactions in the model. If Eq. (7) is not a satisfactory workable empirical model, we could have had recourse to its level form equivalent represented by Eq. (8*). Because of the non-linear functional form with respect to its arguments, this alternative formulation of the theoretical model is not empirically operational in the context of a dynamic econometric model which also includes a spatial dimension. To our knowledge, spatial econometric techniques have not been yet developed to estimate dynamic non-linear specifications. Instead, a vast amount of literature has been addressed to account for spatial interactions in linear models (Anselin and Bera 1998;LeSage and Pace 2009). In particular, Elhorst (2014) describes a set of Dynamic Spatial Panel Model (DSPM) specifications that can be used in the context of stationary panel data. DSPMs are linear and vary according to the dynamics in space and time of the dependent variable, the regressors and the error term. To estimate our empirical model in the context of a DSPM, we formulate a linearized specification, corresponding to a first-order Taylor series approximation of Eq. (8*): where y t = c t , and where the expected signs of the coefficients are: 0 ≥ 0 , 1 ≤ 0 , 2 ≤≥ 0 , and 3 ≤ 0 . This empirical approach enables us to estimate a linear model that integrates spatial panel data interactions with the theoretical benchmark of the previous section.
We construct a panel dataset with Swedish counties as cross-sectional units spanning from 2007 to 2015, described in detail in Sect. 5, where we consider the harvests of three game species as dependent variables: Roe deer, moose and wild boar. As a reference point, we begin by estimating the non-spatial model (Eq. 9) by ordinary least squares (OLS) for each game species. Because ignoring spatial autocorrelation leads to inconsistent OLS estimates and erroneous inference (LeSage and Pace 2009), we extend this baseline model (Eq. 9) to include spatial interaction effects. To this end, we follow the conventional approach in applied spatial econometric studies (Marbuah et al. 2021;Kim et al. 2003) which compare the non-spatial model with different DSPM specifications. The selection of DSPMs addressed in this study is substantiated throughout this section. Because there is a relatively small number of time periods ( N > T ), we set up DSPMs with cross-sectional (county) fixed effects as developed in Elhorst (2014).
To include spatial interaction effects in Eq. (9), we start by formulating a general spatial specification that nests the family set of spatial regression models: where the cross-sectional unit is the county n for time periods t = {2007,…,2015}. The dependent variable y t is a column vector n × 1 , in which n = 20 for the cases of moose and roe deer, whereas n = 14 for wild boar. 8 The scalar is the response of the lagged dependent variable, assumed to be restricted to (−1, 1) to ensure the existence of a long-term stable and convergent equilibrium path (Elhorst 2014). The parameter measures the degree to which the harvest per unit effort is affected by the harvest of neighboring counties. The n × n weighting matrix W characterizes the degree of spatial dependence, and we define it as a first-order contiguity matrix with elements w ij = 1 if i and j are neighboring counties and equal to zero otherwise. As commonly done in the literature, we standardize the weights in each row so that the sum of the row elements equals one. The variable X t−1 is a n × k matrix of (lagged) independent variables, with k = 4 regressors. The four regressors are the lagged hunting effort E t−1 , and the lagged population number of three large carnivores that prey on ungulate game in Sweden: Wolf, brown bear, and lynx. The three lagged carnivore populations correspond to the term L t−1 , which are estimated as separate explanatory variables. The variable is the 4 × 1 vector of partial effects of each of the four regressors. The variable is a 4 × 1 vector of parameters that measures the sensitivity of the harvest to changes in the exogenous variables, that is, the extent to which carnivore numbers and hunting effort in a given county could influence the harvest in contiguous counties. The variable refers to fixed effects across counties, which controls for all omitted variables remaining constant during the nine-year period such as, e.g., landscape features. Conversely, u t represents the time-variant omitted variables such as, e.g., weather conditions, and v t are residuals which are assumed to be normally distributed with zero mean and constant variance 2 . The parameter measures the extent of spatial correlation in the omitted variables.
Due to problems of identification (Manski 1993), we do not estimate the most general spatial model expressed by Eq. (10). Instead, we reduce it to four (nested) spatial specifications: the Spatial Durbin Model (SDM), the Spatial Autoregressive Model (SAR), the Spatial Error Model (SEM) and the Spatial Autoregressive Combined model (SAC or SARAR). The SDM arises by imposing = 0 , and it subsumes the SEM and the SAR models. Namely, the SDM can be reduced to a SAR model by imposing the restriction = 0 , or conversely, it can be reduced to a SEM by imposing = − . On the other hand, the SAC model emerges with a constraint = 0 , which can be reduced to a SAR model by imposing the constrain = 0 , or to a SEM by constraining = 0. Our motivation to compare these four spatial specifications is that it allows us to inspect significant spatial spillovers either in the dependent variable, the independent variables and/or the residuals.
Our data covers the counties where game harvests are different from zero in at least one year during the period 2007-2015, therefore the sample represents the population of counties and not a random sample, in which case the estimation of spatial panels with fixed effects is more appropriate (Beenstock and Felsenstein 2010). We estimate a non-spatial model and four dynamic spatial specifications (SAR, SDM, SEM and SAC) for each of the three game species. The non-spatial model is produced with fixed effects (FE) estimators, while the dynamic spatial models are estimated with a quasi-maximum likelihood (QML) method using a mixed estimator for short time samples. This mixed estimator is based on a bias corrected maximum-likelihood (ML) estimator for the spatial parameter in Wy t , and on an unconditional ML estimator for the other parameters in the model (Elhorst 2010).
In order to know which specification best fits the data, we follow the strategy described in LeSage and Pace (2009). We begin with the SDM and test the restrictions that derive the SAR and SEM. Failing to reject H0: = 0 , with ≠ 0, would favor a SAR specification. On the other hand, failing to reject H0: = − would favor the SEM specification. Because the SAC and SDM are non-nested, the Log-likelihood function value (LL), the AIC (Akaike's Information Criterion) and the BIC (Bayesian Information Criterion) are used to choose the most appropriate model (Belotti et al. 2017;Elhorst 2014).

Marginal Effects
The derivative of the dependent variable with respect to each regressor represents an associated marginal effect, namely the marginal impact of carnivores and/or hunting effort on the game harvests per unit effort. Depending on the dynamic spatial model, it is possible to calculate direct, indirect and total marginal effects for every explanatory variable. Moreover, it is possible to disaggregate these effects into short and long term. However, the SEM and FE estimations only allow for the computation of direct marginal effects (Golgher and Voss 2016).
In the results section, we identify the dynamic spatial-lag model (SAR) as the best fit for roe deer harvests, the dynamic spatial error model (SEM) for wild boar and the non-spatial model (FE) for moose. We therefore restrict the exposition of marginal effects to these models. The detailed procedure to calculate the marginal effects for the three spatial specifications is included in Appendix A.

Addressing Potential Endogeneity with System GMM
The estimated parameters in the dynamic SAR model for roe deer are biased and inconsistent. This is especially true because the spatial lag term ( Wy t ) is correlated with the error term, a classic case of endogeneity. The time-lagged dependent variable ( y t−1 ) on the RHS also induces the so-called Nickell bias (Nickell 1981) since it could be correlated with the county fixed effects. The remaining controls, X t−1 = (wolf t−1 , lynx t−1 , bear t−1 , effort t−1 ) , are potentially endogenous in the model and are likely also correlated with the county fixed effects and the residual term. Purging the fixed effects through first differencing does not entirely solve the identification challenges in the model. This is because the remaining differenced variables, for instance, might still be correlated with the residuals.
These econometric challenges thus render both the pooled OLS and ML estimators the least appropriate for purpose in this study. In order to consistently estimate unbiased model coefficients, we re-estimate the SAR dynamic version of Eq. (10) using an augmented system generalized method of moments (sys-GMM; Blundell and Bond 1998). Kukenova and Monteiro (2009) 9 show that the spatial sys-GMM model outperforms other comparator estimators such as the ML and QML by consistently estimating the spatial lagged coefficient in addition to enhanced efficiency (lower bias) in finite samples and persistent data (see Blundell and Bond 1998;Elhorst 2012). The spatial sys-GMM not only corrects for the endogeneity of the spatially weighted dependent variable, but other covariates considered potentially endogenous. The sys-GMM approach proceeds by estimating the model through first differencing to remove the fixed effect term. Sys-GMM uses internal lagged level and differences of the variables as (valid) instruments (Blundell and Bond 1998;Arellano and Bond 2001). In our roe deer application, we use the second lagged ( roe t−2 ) variable to instrument for roe t−1 and W 2 × roe t for W × roe t . Furthermore, we use the first and second order spatially weighted variables ( WX t−1 and W 2 X t−1 ) for X t−1 = (wolf t−1 , lynx t−1 , bear t−1 , effort t−1 ) following Kelejian and Prucha (1999;2004).
We check the consistency of our adopted estimator by reporting the second order serial correlation test (Arellano and Bond 2001), and the Hansen J overidentification test.

Data
We use data on county level of game harvests, large carnivore populations and number of hunters. For the dependent variables (c t ) , we use hunting bags 10 of roe deer, wild boar and moose. Roe deer and wild boar harvests correspond to the annual hunting bag statistics, collected and administered by the Swedish Hunters' Association (Viltdata 2020). Moose harvest is the annual cull reported by the County Administrative Boards (CAB 2020) for moose management. To run estimations on roe deer and moose harvests, we construct a balanced panel data comprising nine periods T = 2007 − 2015 for N = 20 9 Kukenova and Monteiro (2008) applied this model to study the impact of the differences in environmental regulations on bilateral foreign direct investment (FDI) flows between OECD and developing countries. They find among other things a negative association between environmental stringency and FDI flows after accounting for potential endogeneity and spatial spillovers using the system GMM methodology. 10 A hunting bag is the number of game killed in hunting activities.
Swedish counties ( Figure S1 in the supplementary appendix). 11 For the estimations on wild boar harvest, we only consider N = 14 counties in order to run balanced panel regressions for this species. 12 Hence, the total number of observations available for wild boar is N × T = 126 , while for roe deer and moose, it is N × T = 180 . Note that the dependent variables are defined as harvest per unit effort, therefore harvests (i.e., bags) are divided by the hunting effort variable which will be explained.
For the explanatory variables related to carnivore presence (L t−1 ) , we use information on the estimated populations of wolf, lynx and brown bear. Wolf and lynx inventories were obtained from the reports provided by the Swedish Wildlife Damage Centre (Viltskadecenter 2020). Wolf inventories in Scandinavia are based on snow tracks, droppings and geo-referenced collars from wolf packs and stationary individuals. Lynx inventories relate to the yearly counts of family groups, which are multiplied by an extrapolation factor to determine the population size, see Andrén et al. (2002).  Kindberg and Swenson (2014). From the brown bear inventories, we subtract the licensed bear killings, reported by the National Veterinary Institute (SVA 2020).
We use data on number of hunting licenses per county and year, obtained from the SEPA (Eriksson et al. 2018), as a proxy for hunting effort (E t−1 ) . It is necessary for a hunter to acquire such a hunting license from the SEPA to be eligible to hunt and, thus, one hunting license corresponds to one eligible hunter. The licenses are valid for one year, and the cost for obtaining a license is minor, but it is necessary to have completed a course in wildlife management, and to have proven shooting skills, to qualify for a license.
Due to the large size differences among the Swedish counties, all dependent and explanatory variables were divided by the county area in squared kilometers. County size was retrieved from Statistics Sweden. Table 1 describes all the variables used in this study.

Results
Because all the models include time-lagged covariates (T − 1 = 8) , the regressions are estimated with a total of 112 observations for wild boar harvest, and 160 observations for roe deer and moose harvests. The dependent variable is the contemporaneous harvest per unit effort of the game species at issue c t = h t E t . The standard errors of the spatial models are obtained using Monte Carlo simulations (LeSage and Pace 2009), while the non-spatial models are estimated with conventional standard errors. Results are shown in Table 2. Each column corresponds to a model specification. The FE columns relate to the (non-spatial) fixed effects model, while the rest of the columns correspond to the dynamic spatial specifications addressed in this study (SAR, SEM, SDM and SAC). Columns (1)-(5) display the results for roe deer harvests per unit effort, columns (6)-(10) for boar harvests per unit effort, and columns (11)-(15) for moose harvests per unit effort. To validate our empirical estimations, we assessed the robustness of the results in Table 2 to alternative weighting matrices. The robustness checks have been included in Appendix C.

Roe Deer Harvests
The high and positive significance of the time-lagged roe harvest coefficient ( ) in all specifications provides evidence of dynamic effects. The positive sign of the coefficient is intuitively reasonable as it can be associated to natural variability in the game populations. The high significance suggests that game harvests are time-dependent, and as a result, the roe deer harvest in the preceding year can be used to predict current harvest values. From all regressions, lynx is the only regressor to have a direct negative impact on roe deer harvest at the 5-10% level of significance. The lynx point estimates reveal an upward bias (in absolute terms) of the lynx coefficients when ignoring the spatial autocorrelation. This overestimation can go from 7 to 46% based on a comparison of the lynx coefficients of the spatial specifications with that of the non-spatial linear model. The significant spatial autoregressive parameter ( ) for the SAR, SDM and SAC models show evidence of spatial autocorrelation occurring in the dependent variable amid cross-sectional units. Namely, there are roe deer harvest spillovers into other counties that are statistically significant, and therefore justify the use of spatial econometric modelling. The positive sign indicates that game harvests may increase in neighboring counties as ungulates migrate trying to reduce  predation risk (Hebblewhite and Merrill 2007;Cromsigt et al. 2013). These positive spatial interactions in harvests are also consistent with the unregulated harvest levels, and the fact that hunters can choose harvest levels based on their preferences. The existence of spatial interactions is corroborated with the results for the spatial parameter of the residuals ( ) , which is positive and significant for both the SEM and SAC specifications. The specification tests for SDM fail to reject that = 0( 2 = 0.78) , hence favouring the SAR model. Similarly, the SEM is favored over the SDM provided that the null hypothesis = − cannot be rejected ( 2 = 0.08) . However, the SAR is slightly but systematically superior to the SEM when comparing the log-likelihood (LL) function values, the AIC and *** p < 0.01, ** p < 0.05, * p < 0.1. Standard errors in parentheses. We do not include the spatially lagged explanatory variables (WX t−1 ) from the SDM, because all coefficients ( ) are not statistically different from zero. The results for the spatially lagged carnivore regressors have been omitted from the SDM column given that none are statistically different from zero the BIC. Similarly, the SAC model is preferred to the SDM given the higher LL and lower AIC and BIC. These information criteria yield rather similar results for the SAC and SAR specifications. The SAR and SAC have the same number of explanatory variables, hence the R 2 is appropriate to compare the fit of the two models. The SAR model outperforms the SAC when accounting for the variation between panel units (58.27%), within panel units (15.91%) and overall (49.9%). Provided the statistical and behavioral reasons stated above, we choose the SAR specification to best describe roe deer harvest data.
To calculate the partial effects, we estimate the SAR model using a spatial dynamic sys-GMM estimator. The use of the sys-GMM method is necessary to identify the SAR coefficients, which are a combination of direct and indirect effects (StataCorp 2019). We proceed with the estimation using different weighting matrices (Table 3).
We choose the model capturing the largest variation (i.e., with highest adjusted R 2 ), which corresponds to the sys-GMM model using inverse distance spatial weight matrix (model 2). The selected model has the least number of instruments in comparison to Table 3 Estimation of the spatial dynamic sys-GMM model using different weighting matrices Dependent variable: Roe deer harvest per unit effort. Robust standard errors clustered at the county-level in parentheses. *** p < 0.01 (1%), ** p < 0.05 (5%), * p < 0.1 (10%). Results in columns (3)-(5) are based on k-nearest neighbour(s) spatial weight matrix, where k = 1,2 and 3 neighbours respectively. AR(2) implies residual autocorrelation of the 2 nd order in each model whereas Hansen J and Difference-in-Hansen test statistics respectively correspond to testing the null hypotheses of validity and exogeneity of the instrumental variables. Computation of marginal effects is based on our preferred model (2) Table A1 (for the SAR model) using the parameters ( , , lynx ) obtained in column (2) of Table 3. The partial impacts show evidence of significant indirect effects in the short and long run. That is, lynx species seem to create both short-and long-term spillovers on roe deer harvests. The point estimate of the non-spatial model ( −16.98 * * ; column (1) of Table 2) is four times larger than the short-term marginal impact of roe deer ( −4.39 * * ; column (2) of Table 3). Namely, the non-spatial model overestimates the marginal effect compared with that of the spatial specification. Thus, if policy makers take the impact on roe deer hunting into account, and strive towards an economically optimal lynx population, the lynx population level will be lower than optimal if the decisions is based on results from non-spatial studies.

Wild Boar Harvests
The results for wild boar yield highly significant dynamic effects in all specifications, that is, the time-lagged boar harvest coefficient ( ) is always significant and positive. Similarly, the highly significant parameters of spatial autocorrelation of the dependent variable ( ) and of the disturbances ( ) indicate the presence of spatial influences in the harvest of this species. The highly significant coefficient can be explained by relevant omitted variables contained in the residuals, e.g., landscape characteristics and habitat conditions influencing the growth rate of boar populations . Our findings indicate that hunting effort becomes significant only when the spatial structure is considered. The lagged hunting effort is significant at 5% level in all dynamic spatial models ( Ef fort ) , as well as the spatial lag of effort in neighboring panel units ( Ef fort ) . The former yields negative signs suggesting that the harvest decreases directly due to the hunters in the county. The latter positive outcome (1.027 *** ; column (8) of Table 2) indicates that the harvest can also increase indirectly via incoming migration of wild boar caused by hunting efforts in neighboring counties. None of the large carnivore species has a statistically significant impact on boar harvest either directly or indirectly, which is consistent with expectations.
Proceeding in a similar fashion as with roe deer harvests in Sect. 6.1, the performance tests yield the SEM as the best fit for the boar harvest data. The selection of this model implies the calculation of the direct marginal effect in the short and long terms (Belotti et al. 2017).

Moose Harvests
As opposed to roe deer and wild boar estimations, the results for moose are non-stable when considering the dynamic component in the spatial models. This is because the stability condition + < 1 is not fulfilled when running the regressions (Elhorst 2014). 13 As a result, our approach is to analyze the moose harvests assuming = 0 in Eq. (10), thus consider the spatial dimension without the dynamic component. Table A2 in Appendix A shows the results for moose using the spatial dynamic sys-GMM. The results reinforce our choice of the spatial FE model for this species, because no evidence is found of statistically significant spatial spillovers in the dependent variable or in the residuals for any of the spatial specifications. This translates into relatively better performance indicators for the non-spatial model with the LL measure as the only exception. We also estimated a bias-corrected non-spatial dynamic panel fixed effects model for each species (Table A3, Appendix A). The interested reader can refer to Everaert and Pozzi (2007), De Vos et al. (2015), and Kiviet (1995 for further details. The lack of dynamic and spatial effects in harvests could be related with the use of harvesting quotas, as the administrative decision process may not be flexible enough to fully adjust quotas to changes in moose abundance due to migration. The lagged wolf covariate yields a negatively significant effect at 1% level in all cases, which is intuitively consistent given that wolves are effective preying on large ungulates (Persson 1996). The lagged bear parameter is negatively significant at 10% level in the non-spatial specification. The brown bear is significant at the 5% level in the spatial estimations, and the SDM specification yields some positive indirect effects from bears of neighboring counties (at 10% level).

Calculation of Marginal Effects
For policy relevance and intuitive purposes, we express the marginal effects as that of one additional unit of carnivore or hunter on units of game harvested. We use the formulas in Table A1 and the parameters in Table A4 to determine the marginal effects for each game species. By holding constant the effort variable and the county area at their respective average values E = 13, 198.1;Countyarea = 20, 390 Table 2, we can compute the average impact of an additional unit of the regressor (i.e., carnivore or hunter) on units of game bags, for details on the calculations see Appendix B. Table 4 shows the resulting marginal effects.

Marginal Implicit Prices of Carnivores and Hunting Effort
Based on the marginal effects calculated in Table 4, we calculate the marginal implicit price of an additional carnivore (or hunting effort) in terms of the impact on game harvests.  Mattsson et al. (2008) estimates that the gross hunting value of roe deer in the hunting season 2005/6 was about 712 million Swedish kroner (SEK), or equivalently, 78 million US dollars (USD). That year the roe deer harvest was 129,000 (www. viltd ata. se). The average hunting value per roe deer can then be estimated to SEK 5,519 (USD 610). We use this value to estimate the implicit price of lynx in terms of its impact on roe deer hunting. Mensah and Elofsson (2017) estimate that the average marginal willingness-to-pay for an extra unit of wild boar harvest is 2,709 SEK (299 USD). Provided our result and the latter study, an additional hunter (i.e., hunting license) reduces boar hunting value by 948 SEK (104 USD) and by 1,354 SEK (149 USD) in the short and long run, respectively.
According to Boman et al (2011), the average willingness-to-pay per harvested moose is 9,207 SEK (1,018 USD). Multiplying this with our estimated marginal effect implies that one additional wolf decreases moose hunting value by 139,946 SEK (15,480 USD). For brown bear, the same procedure yields 38,116 SEK (4,216 USD). The resulting values are presented in Table 5.

Discussion and Conclusion
We apply a novel approach to wildlife management and conservation, which captures the time dynamics of population development as well as spatial spillovers due to species' and human mobility. Building on a bioeconomic model of a predator-prey system, with the prey being harvested, we apply spatial econometric techniques to quantify the impact of human hunting, and predation by three large carnivores, on the harvests of three major game species, as well as spatial spillovers in harvests. Our study contributes to the literature by quantifying the magnitude of spatial spillovers for a harvested prey species, and by comparing those across three game species for which the harvesting regulation differ.
Our results reveal the presence of dynamic effects and spatial spillovers in the harvests of roe deer, where harvest levels are unregulated, and wild boar, where harvesting is encouraged due to the damage on agricultural crops. In contrast, results show no strong evidence of spatial spillovers in the moose harvests, consistent with the fact that moose harvesting is regulated at the regional level. If harvest opportunities decline within a Moose Management Area, hunters cannot travel to another Moose Management Area and increase the harvest there above the quota level. Intense predation within a county could, hypothetically, create moose populations to shift to nearby regions. This could allow for larger quotas in the neighbor regions. However, ecological research suggests wolf presence does not affect moose' choice of habitat and location (Sand et al. 2019). Also, it might be that the administrative process to reach a decision on harvest quotas responds slowly to moose population changes caused by migration. Both of those could contribute to low levels of spillovers.
After accounting for the dynamic spatial regimes, lynx predation is found to have a significant negative influence on roe deer harvests in the same county in the short-term (a direct effect) as well as in the neighboring counties (an indirect effect). This is consistent  with the fact that the Eurasian lynx mainly preys on roe deer and therefore influences the roe population dynamics (Andrén and Liberg 2015). It is also consistent with lynx migrating to find prey. Our estimates of the impact of lynx on roe deer hunting values could be compared to hedonic pricing estimates of the impact of lynx on hunting lease in Mensah et al. (2019) and Lozano et al. (2020). The former study concludes that the marginal implicit price of an additional lynx territory is SEK 1.51 million (USD 180 thousand), whereas the latter estimates the same effect to be SEK 1.31 million (USD 145 thousand). Using our results above and applying a conversion factor of six (Andrén et al. 2002) to account for the average number of lynx individuals in a territory, the total effect of an additional lynx territory in the short run would amount to SEK 0.16 million (USD 18 thousand), which is smaller than the earlier estimates. One potential explanation is that the present study does not account for impact of lynx on small game. However, another possibility is that the production function approach might yield lower costs estimates than the hedonic pricing method.
For wild boar harvests, there is significant evidence of spatial interactions in unobserved characteristics, suggesting the presence of spatial patterns not captured by hunting pressure or by carnivore abundance. Higher human hunting effort reduces boar harvests in the same county, indicating that deliberate efforts to reduce the considerable damages to agricultural crops by this species (Gren et al. 2020) have had some success. Such efforts can be driven by communication from national and regional wildlife management agencies, as well as local pressure from landowners on hunters that lease hunting land.
Our results show that an increase in wolf numbers limits locally the possible harvests of moose, which is in line with results in Boman et al. (2003) and Lozano et al. (2020). The latter estimates the marginal implicit price of a wolf territory, in terms of the impact on hunting lease prices, to SEK 1.5 million (USD 166 thousand). Based on Svensson et al. (2019), the ratio of wolf individuals and wolf territories in Sweden is about 5.8. Hence, the direct effect in this study corresponds to an implicit marginal price of a wolf territory in terms of the impact on moose harvest equal to SEK 0.81 million (USD 89 thousand). This is about half of the cost estimated by Lozano et al. (2020). A lower effect should be expected in the present study, as wolves are large and can also prey on other species (Persson 1996). In spite of human harvest being the main limiting and regulating factor of moose population in Sweden (Ericsson 1999), our results do not provide statistically significant estimates for hunting effort which is presumably due to the harvest quotas which entail that hunters' own preferences have a very modest impact on moose populations.
A limitation of this study is the large size of the cross-sectional units, which is due to the fact that harvest data useful for the analysis is only available at the county level. This implies that the presence of spatial effects can be overlooked compared with more disaggregated data. A smaller spatial aggregation can potentially produce more accurate estimations, but such data are currently not available.
Knowledge on the marginal implicit price of large carnivores, in terms of their impact on hunting, are important for policy because together with information on the marginal benefits, and on other costs, e.g. for depredation on livestock, they can inform on optimal population sizes. Existing environmental valuation studies applied to large carnivores in the Nordic countries using stated preferences methods do not inform on marginal benefits, but on the average willingness-to-pay for conservation. In addition, there are economic values associated with brown bear license hunting (Lozano et al. 2020). Yet, it is evident from the literature that the benefits substantially exceed costs (cf., e.g., Håkansson et al. 2011;Widman and Elofsson 2018). Our results reveal that the marginal implicit prices of carnivores can be substantially overestimated if non-spatial models are applied (such as in, e.g., Bostedt and Grahn (2008) and Widman and Elofsson (2018)). In particular, our study shows that larger populations of lynx will be economically optimal when spatial spillovers are taken into account. Knowledge on the economic impact of large carnivores on hunting values is also important to avoid exaggerations in the debate, and conflict among stakeholders regarding wildlife management.