Abstract
The random utility model (RUM) is commonly used in the land-use and fishery economics literature. This research investigates the affect that spatial heterogeneity and spatial autocorrelation have within the RUM framework using alternative specifications of the multinomial logit, multinomial probit, and spatial multinomial probit models. Using data on the spatial decisions of fishermen, the results illustrate that ignoring spatial heterogeneity in the unobservable portion of the RUM dramatically affects model performance and welfare estimates. Furthermore, accounting for spatial autocorrelation in addition to spatial heterogeneity increases the performance of the RUM. (JEL C35, Q22)
I. Introduction
Within the environmental economics literature, McFadden’s (1974, 1978) random utility model (RUM) has been extensively used when analyzing microlevel discrete choices. Although a number of settings do not explicitly depend upon spatial choice alternatives, there are a number of applications where space is an important element of the choice. Most notable are land-use models (e.g., McMillen 1989; Bockstael 1996) and the location choice models used in the fisheries literature (e.g., Eales and Wilen 1986; Holland and Sutinen 2000). The standard convention is to assume a computationally tractable structure for the unobservables in the econometric model. However, when one further investigates these modeling assumptions it is clear that they are untenable in a spatial discrete choice model. This research investigates these assumptions in detail and proposes an empirical model that accounts for spatial heterogeneity and autocorrelation in the unobserved portion of utility.
The earliest utilization of a spatial discrete choice model to investigate locational choice decisions in fisheries was conducted by Eales and Wilen (1986) in their study of the pink shrimp fishery.1 Since that time these models have been used to investigate locational choice decisions in the New England groundfish fishery (Holland and Sutinen 1999, 2000), the welfare implications of spatial closures (Curtis and Hicks 2000; Dalton and Ralston 2004; Hicks, Kirkley, and Strand 2004; Smith and Wilen 2004), heterogeneous risk preferences (Mistiaen and Strand 2000), and state dependence (Smith 2005), as well as to predict the behavioral responses of fishermen to spatial regulations (Hicks and Schnier 2010; Smith 2004; Smith and Wilen 2003), to cite a few.
In the fisheries literature most of the models characterize decisions regarding where to fish at a given point in time. The choice set possessed by the decision agent is explicitly spatial. Therefore, it is plausible that the unobserved portion of site-specific utility in a RUM exhibits spatial heterogeneity as well as autocorrelation. Spatial heterogeneity in models of fishing behavior may arise due to unobserved benthic habitat characteristics, stock densities, and nutrient-rich waters, which are location specific and known to the decision agent but not to the researcher (Hicks and Schnier 2010; Smith 2010). Spatial correlation in the unobserved portion of utility in a RUM can arise from two primary sources: (1) spatially correlated omitted variables and (2) asymmetries in the definition of “space” assumed by the researcher versus the “true” definition of space in which the decision-agents operate (Bell and Irwin 2002). Spatially correlated omitted variables in models of fishing behavior may arise due to correlations in the site-specific characteristics (e.g., spatial gradients in the benthic characteristics). Spatial correlation due to space asymmetries may arise when the spatial grids a researcher uses to compartmentalize effort do not match up with the spatial grids fishermen use in their decision of where to fish. In fact, most spatial discrete choice modeling uses an arbitrarily defined definition of “space” that is rarely perfectly correlated with the spatial definition utilized by the decision agent. Therefore, it is plausible to assume that many spatial discrete choice models are subject to the presence of spatial autocorrelation. The presence of spatial autocorrelation in discrete choice models will yield inconsistent and inefficient parameter estimates.
The general results from this research indicate that a failure to account for spatial heterogeneity in a discrete choice model compromises the performance of one’s estimator. Accounting for spatial heterogeneity can be accomplished using either a multinomial logit (MNL) model with alternative specific constants (ASCs) or with a multinomial probit (MNP) model that incorporates site-specific Cholesky factors. Furthermore, accounting for spatial autocorrelation within the MNP model (estimating a spatial multinomial probit [SMNP]) will further increase the performance. The following section outlines the empirical models estimated that account for varying degrees of spatial heterogeneity and autocorrelation in unobserved site-specific utility.
II. Empirical Model
Consider a fisherman f deciding among i = 1,..., N different locations in which to fish. Assuming a linear-in-parameters utility specification, the utility they derive from fishing in location i can be expressed as
[1]The observation matrix Xif contains information that varies across fishermen and locations, β is a preference parameter assumed to be constant across all fishermen in the population, and εif is the portion of site-specific utility that is assumed to be known by the decision agent (fisherman) but unknown to the researcher.2 The distributional properties of the vector of unobserved utility ε = {εif,...,εNf} is the primary focus of this research, where f (ε) represents their joint density. The probability that a decision agent (fisherman) will choose the spatial alternative i is represented as
[2]Using the joint density of ε the probability can be expressed as (Train 2003)
[3]where I( · ) is binary indicator function determining whether or not the spatial alternative was selected. The estimation of the behavioral parameters in equation [3], β, requires the estimation of a multidimensional integral over f (ε), which depending on its specification yields different types of discrete choice models. However, there are only a few types of distributions over which the model has a closed form solution. This is precisely why many researchers assume that the errors are independent and identically distributed (iid) Type I extreme value (Gumbel distribution). This assumption facilitates the estimation of the choice probabilities via a logit (or MNL when N > 2) specification. Given this specification the choice probability can be represented as (Train 2003)
[4]where σ is the variance scale parameter that is not identified in the empirical model. Therefore, the researcher estimates β/σ, where the implied variance of the error structure is σ2(π2/6) for all alternatives in the choice set (Train 2003). This model has become the workhorse for many empirical applications within the spatial choice literature (Eales and Wilen 1986; Holland and Sutinen 2000; Hicks, Kirkley, and Strand 2004). Although alternative specifications of the error density f (ε) yield a more robust model of spatial correlation, they come at a higher computational cost because estimation requires the simulation of the error density f(ε). Three alternative specifications of the unobserved utility will be investigated in this research: an MNL with ASCs, an MNP, and an SMNP. The MNL with ASCs utilizes site-specific fixed effects to control for time-invariant spatial heterogeneity. The MNP model controls for spatial heterogeneity using site-specific error variances, and the SMNP model controls for spatial heterogeneity and spatial autocorrelation.
Multinomial Logit with Alternative Specific Constants
Although the MNL model does not explicitly account for spatial heterogeneity and autocorrelation within the error structure, it is possible to account for spatial heterogeneity using ASCs. Incorporating ASCs expands equation [1] to
[5]where δi is the site-specific ASC (fixed effects) used to capture latent spatial heterogeneity. ASCs have been used in the fisheries literature to control for unobserved site-specific utility measures such as the biological characteristics of the fishery (Hicks and Schnier 2010) and to investigate congestion/agglomeration effects in commercial (Hicks, Horrace, and Schnier 2009) and recreational fisheries (Timmins and Murdock 2007). Estimation of equation [5] is achieved using the contraction mapping proposed by Berry (1994) and utilized by Berry, Levinson, and Pakes (1995) in their research on product differentiation. This model will also be estimated in our comparison of the different empirical methods to investigate its relative performance because it is an easy way to account for spatial heterogeneity using the MNL framework.
Multinomial Probit
Our illustration of the MNP model closely follows that of Train (2003) and Bolduc (1999). The MNP model assumes that the unobserved portion of site-specific utility ε ~ N(0,Σ), where
[6]Given this specification, the spatial correlation in the unobserved site-specific utility is captured by the covariance parameters σj∀i ≠ j, and spatial heterogeneity is captured by 
. However, estimation of Σ is complicated by two factors: (1) the multidimensional integral captured by f (ε), and (2) the fact that as the number of spatial alternatives increases it may be infeasible to estimate all the covariance parameters σij within the model. The former complication can be addressed using the GHK simulator (Geweke 1989, 1991; Hajivassiliou, McFadden, and Ruud 1996; Hajivassiliou and McFadden 1998; Keane 1994). The latter complication is addressed in the specification of the SMNP model.
Estimation first proceeds by normalizing the model for scale and level. This is achieved by expressing each site-specific utility function as the difference between itself and the normalized site-specific utility level (location one in the empirical application), and by setting the variance of the second site-specific error, 
equal to one (Train 2003).3 The transformation redefines the site-specific utility functions as
[7]where ~ N(0,
) and 
is derived from the underlying variance-covariance matrix Σ (Train 2003). The Cholesky matrix for 
can be defined as
[8]The variance-covariance matrix 
can be defined as 
. The extra row and column of zeros added to the matrix are utilized to facilitate estimation (see Train 2003 for more discussion). The parameters that will be estimated are the Cholesky factors κi,j contained within Cε, in addition to the behavioral parameters, β. Given the large number of spatial alternatives in the empirical application it is assumed that κij = 0 ∀ i ≠ j and the SMNP will be used to capture spatial correlation across sites.
Estimation of the choice probabilities via the GHK simulator proceeds by first differencing the transformed utilities, 
, relative to the observed choice. Therefore, the probability that a site is selected, Pif, becomes an integral over the nonpositive orthant (Bolduc 1999; Train 2003).
The simulated probabilities are estimated K times to construct a simulated probability,
[9]which is log-transformed and then maximized to obtain estimates of the behavioral parameter β and the elements of Cε .4 Within the empirical analysis utilized in this research K = 100; alternative simulation levels were experimented with, but when K exceeded 50 the parameter estimates stabilized. Although a higher level of K would reduce the estimation bias, it also comes at a price of higher computational time; therefore, 100 was selected as a reasonable level of simulation.5 Having outlined the MNP and estimation strategy, the following section builds on the MNP to construct the SMNP model.
Spatial Multinomial Probit
In the case that the number of spatial alternatives is large, as will be the case in the empirical example, it may become intractable to estimate the complete set of off-diagonal elements contained in Cε that capture the spatial correlation in the unobserved portion of site-specific utility. In the presence of spatial correlation in the unobserved portion of site-specific utility, equation [1] becomes
[10]where W is a spatial weighting matrix and μ ~ N(0,Σμ). Equation [10] can be also be expressed as
[11]where I is the identity matrix of dimensions N × N. Defining 
and 
Equation [11] can be written as
[12]This specification is similar to the model developed by Bolduc (1999) where the W matrix was used to capture correlations in transportation alternatives. Given this specification, the implied distribution of εif is εif ∼ N(0,(I − ρW)-1 Σμ(I − ρW)-1). In the case that there is no spatial correlation in error structure, the distribution of εif becomes εif ∼ N(0,Σμ), which is identical to the distributional assumption in the MNP because the variance and covariance parameters contained in Σμ are identical to those contained in Σ. Estimation of the model proceeds in a similar fashion to the methods outlined in the MNP section. The utility functions, 
, are normalized for scale and level by differencing each utility by the first site alternative, and then the GHK simulator is used to simulate the probabilities of a site being visited. Denoting 
as the variance-covariance matrix following model normalization, the primary difference between the model outlined in the MNP section and this section is that the Cholesky factor matrix for the SMNP contains only the diagonal elements of Cε, denoted Cμ. Using Cμ the variance-covariance matrix for the model depicted in equation [10] can be expressed as 
, where the Cholesky factor matrix is
[13]Therefore, following estimation of the elements contained in Cμ, it is possible to estimate the spatially correlated elements of unobserved site-specific utility by transforming the estimated variance-covariance matrix by (I − ρW)-1 because 
. This substantially reduces the number of parameters that must be estimated in the model when compared to the full variance-covariance matrix captured by Cε in the MNP. However, it is possible that forcing all spatial correlation to be accounted by the weighting matrix and a single parameter, ρ, is restrictive versus the more flexible specification containing the σij's
III. Empirical Application
The empirical application selected for this analysis uses data from the Atka mackerel fishery in the Bering Sea and Aleutian Islands (BSAI). Atka mackerel is a semipelagic schooling species that is regionally concentrated within the BSAI near the islands and seamounts contained within the island chain and is rarely found in the open ocean (Lowe et al. 2006). Although Atka mackerel is also found in the Gulf of Alaska, conventional wisdom is that these two stocks are spatially distinct (Lowe et al. 2006). Commercial fishing for the Atka mackerel began in the 1970s and was primarily conducted by Russian, Japanese, and Korean vessels up until the inception the Fishery Conservation and Management Act of 1976, which was later renamed the Magnuson Stevens Fishery Conservation and Management Act (MSFCMA). The MSFCMA shifted the fishery to a joint venture fishery up until 1989 when the fishery effectively became a domestic fishery (Lowe et al. 2006).
Atka mackerel is an integral food source for Steller sea lions (Merrick, Chumbley, and Byrd 1997; Guénette et al. 2006), a species protected under the Endangered Species Act, and is spatially regulated to limit fishing activity near the Steller sea lion rookeries (breeding grounds). The North Pacific Fishery Management Council (NPFMC) divides the BSAI into three regional spatial zones defined as the western, central, and eastern zones. These zones are contained in the National Marine Fisheries Service (NMFS) statistical management zones 543 (western), 542 (central), and 541 (eastern), and each zone contains its own regional total allowable catch (TAC) (Lowe et al. 2006). Harvesting within the Atka mackerel fishery is also temporally divided into an A season and a B season. The
A season typically runs from January 1 to April 15, and the B season runs from September 1 to November 1.
The NMFS statistical management zones 542 and 543 (central and western regions) are also further divided into those regions that are critical habitat for the Steller sea lions, referred to as harvest limitation areas (HLAs), which are defined as the waters west of 178 degrees longitude and within 20 nautical miles of an island listed as a Steller sea lion rookery (CFR §§ 679.20). Figure 1 provides a graphical illustration of the area studied, as well as the HLA zones contained within the BSAI. The HLA regions are some of the most productive regions within the Atka mackerel fishery (Hicks and Schnier 2010; Horrace and Schnier 2009). To allow fishermen access into these regions the regional TACs are divided up so that 40% is caught outside of the HLA region and 60% is caught within the HLA. However, to ensure that harvesting is not temporally concentrated in any one fishing region, fishermen are randomly assigned into one of two “platoons” each season. A platoon consists of a set of vessels that have access rights to one of the HLA regions contained in either 543 or 542, with the other platoon assigned to the other region. The regional TACs are also evenly divided among the two platoons, and the platoons switch locations after both have caught all of their regional allocations. Platoon selection determines which spatial locations a fishermen can consider at a given point in time (Hicks and Schnier 2010), and efforts will be made to control for this in the empirical analysis.6
National Marine Fisheries Service Management Zones (Three-Digit Site Designations) and HLA Zones (Light-Shaded Circles around the Islands within the Aleutian Island Chain)
The data set used in the analysis was obtained from the NMFS Observer Program Database. The data set contains information at the haul level (gear deployment) and records the latitude and longitude coordinates of the haul as well as its catch composition and other spatial information such as the management zone in which the haul took place and whether or not a marine mammal interaction took place. Using this information each haul is assigned to one of the management zones with zones 542 and 543 divided into non-HLA (5420 and 5430) and HLA zones (5421 and 5431), as well as to the Alaska Department of Fish and Game (ADF&G) statistical reporting zone. The AFD&G statistical reporting zones divide the BSAI into 1 degree latitude by ½ degree longitude squares and are a smaller spatial resolution than management zones 541, 542, and 543.7 The ADF&G zones were also further divided into non-HLA and HLA zones when the zone contained HLA and non-HLA regions.
The larger management zones effectively divide the BSAI into 6 spatial zones,8 and the ADF&G zones divide effort into 24 distinct spatial zones.9 In addition, to the 6-site and 24-site model, a 12-site model was constructed by aggregating subsets of the 24 sites. Within the 12-site model, sites were longitudinally aggregated except for when a site was contained within the HLA. Table 1 lists the descriptive statistics for 4,287 observations on the 11 vessels that participated within the Atka mackerel fishery from 2002 to 2006 over 100 unique “trips.” Trips were constructed by identifying substantial time gaps between fishing activity that indicated that a vessel had returned to port. This was conducted because the precise start date of each trip is not known with certainty. To reduce the estimation bias resulting for this informational gap, the first haul from each of the trips was removed. Therefore, the empirical models estimated are conditioned on the vessels being on the fishing grounds versus leaving port and then selecting where to fish.10 The effects of this assumption will be discussed in more detail within the results section.
Descriptive Statistics
The average trip level revenues were over a half of a million dollars, and an average trip consisted of nearly 43 hauls.11 In the 6-site model, fishermen visited an average four sites, whereas in the 12-site and 24-site models they visited approximately five and six sites, respectively. Therefore, there exists a moderate degree of spatial movement at the trip level within the Atka mackerel fishery. Using these three definitions of space, the MNL, MNL with ASCs, MNP, and SMNP models will be estimated and compared. Given that one primary focus of this research is the out-of-sample predictive accuracy of the different models, the empirical models will be estimated using data from 2002 to 2005. Data for 2006 will be used to test the out-of-sample predictive accuracy of each model. Furthermore, to facilitate the moving average calculation used to estimate expected revenues, we used a “burn in” period of 30 days. This alteration to the data set, combined with the removal of the first haul on each trip, reduces the sample to 2,946 observations for 2002–2005 and 956 for 2006. In addition to comparing the within-sample and out-of-sample behavioral predictions, we will estimate the compensating variation resulting from a hypothetical site closure (site 5431) using the data from 2002 to 2005. This analysis is conducted because welfare analysis is often conducted within the fisheries literature (Curtis and Hicks 2000; Hicks, Kirkley, and Strand 2004), and comparing the models using this metric may have important policy implications. Having outlined the data being used, the following section discusses the specification of the utility function used in the analysis as well as the data contained in the model.
Utility Specification
The independent variables contained in the utility function discussed above closely follow the empirical specification used by Hicks and Schnier (2010) and are
[14]Disti|k is the distance that vessel expects to travel to visit site i conditional on being in site k expressed in kilometers, with site-specific coordinates calculated as the average latitude and longitude of effort conducted within the site. Rvni is the expected site-specific revenues, which is calculated as a one-year moving average of the revenues observed within the site.12 Stelleri is the expected number of Steller sea lion interactions that is measured as a one-year moving count of the number of Steller sea lion interactions that occur within the site. The final variable, MissDumi, is a binary variable that takes a value of one if a given site was not visited in the previous year and zero otherwise.
Given that some sites are not available to the fishermen at given points in time due to the spatial lottery system used to define platoons, it is important to make sure that site-specific utility accounts for this restriction. This was achieved by removing those sites that are unavailable to fishermen at a given point in time from their spatial choice set. Furthermore, the Cholesky factor matrix was altered by removing the rows and columns that represent unavailable sites based on the HLA lottery system. Therefore, only those sites that are contained in a fisherman’s feasible set of sites to visit are given empirical weight in the analysis.
Within the SMNP models the spatial weighting matrix used defines each element of W, defined as ωij, to be the inverse of the distance traveled from the mean reported latitude and longitude of site i to site j squared,13
[15]Alternative specifications of the weighting matrix were utilized and each generated a similar profile of parameter estimates and Cholesky factors.14 However, the inverse distance squared function consistently generated a higher log-likelihood value than the alternative specifications and was therefore selected.15 In addition, the weighting matrix was row-standardized and the spatial autocorrelation parameter, ρ, was restricted to ρ ∈ ( − 1,1). Spatial weighting matrices were not altered to control for whether or not a site was “open” to a vessel at a given point in time, because the expected utility they derive from visiting a site may be spatially correlated with the expected utility they would derive if a closed site became open.16
IV. Results
Table 2 illustrates the modeling assumptions and structural differences across the six empirical models that are estimated for each of the three spatial resolutions: the MNL (MNL), MNL with ASCs (MNLH), MNP with a constant variance (MNP), MNP with spatially varying variance parameters (MNPH), SMNP with constant variance (SMNP), and SMNP with spatially varying variance parameters (SMNPH). The models with spatially varying variance parameters (MNPH and SMNPH) and spatially varying constants (MNLH) control for site-specific spatial heterogeneity, hence the “H” at the end of the their acronym. Spatial autocorrelation indicates whether or not the spatial structure of the unobservables, see equation [10], was accounted for in the estimation routine. This is captured by the SMNP and SMNPH models and represented by an “S” at the beginning of the acronym. The results are reported in Tables 3 through 5. Table 3 contains the results for the 6-site choice model consisting of the three NMFS management zones with regions 542 and 543 further subdivided into non-HLA (5420 and 5430) and HLA zones (5421 and 5431), and region 5000 used to define all other spatial activity. Table 4 contains the results for the 12-site choice model that further divides the NMFS management zones, and Table 5 contains the results when the spatial resolution of the sites is modeled at the ADF&G site resolution, consisting of 24 sites.
Comparison of the Six Empirical Models
Parameter Estimates for the Six-Site Location Choice Model
Parameter Estimates for the 12-Site Location Model
Parameter Estimates for the 24-Site Location Model
Six-Site Model Results
The six-site model results illustrate that the expected distance traveled has a negative and statistically significant effect on the probability that a site is visited across all six models estimated. The expected revenues within a site are not statistically significant for three of the models estimated, but they are statistically significant and positive for MNPH, and statistically significant and negative for SMNP and SMNPH. The expected number of Stellar sea lion interactions is positive and statistically significant for all of the models except MNLH, and the missing dummy variable, which indicates if vessels have not fished in a site within the past month, is positive and statistically significant. Although the coefficient on distance makes economic sense, the coefficients on revenues and the missing dummy variables warrant some additional investigation.
The coefficients on revenues as well as that on the missing dummy variable are a construct of the empirical specification. In order to minimize the estimation bias resulting from not knowing the true start date for each trip, we removed the first haul on the constructed trip. This assumption places the fishermen in one of the six sites and then estimates the probability that they will either remain in the site or leave for another site. The revenue coefficients indicate that fishermen are very reluctant to leave one of the six macroregions as the spatial revenue gradients are often time insufficient to motivate a move; in essence fishermen as very “sluggish” in their spatial responses when we use this very large spatial resolution. However, fishermen do seasonally switch locations, and the coefficient on the missing dummy variable is capturing these seasonal changes.17 As will become more evident shortly, these coefficients behave more consistently with our a priori expectations when we expand the spatial resolution. This said, these models can be used to statistically predict behavior because they assign relative statistical weights to the information present. However, they cannot be used to investigate welfare effects.
Analyzing the Cholesky factors estimated for the six-site choice model (Table 3) illustrates that the variance in the unobservable portion of utility in all sites is less than that observed in Site 2 (5410). This indicates that Site 5410 possesses the highest level of random variation in the site-specific utility. Furthermore, comparing the implied variance-covariance matrix for the random utility model using SMNP and SMNPH, 
, indicates that they are very similar to one another. This is due to the fact that the constant variance parameter estimated in SMNP is similar to the average of the variance parameters estimated in SMNPH. More importantly the utilization of the spatial autocorrelation model expands the model to account for spatial heterogeneity and spatial autocorrelation by estimating an implied covariance across all the sites in the model. Lastly, the results indicate that spatial autocorrelation is present, as ρ is statistically significant in both the SMNP and SMNPH models.
Twelve-Site Model Results
The 12-site model results accord with our a priori expectations. The coefficients on distance and the missing dummy variable are all statistically significant and negative, and the coefficients on expected revenues are all positive and statistically significant. As was the case with the six-site model, the first haul of each trip was also deleted within this data set. The difference in the coefficients between the 6-site model and the 12-site model indicate that at the 12-site spatial resolution spatial mobility becomes more evident and that revenue gradients are sufficient to generate movement within the fishery. This indicates that the “sluggish” response discovered at the six-site resolution dissipates at the slightly larger spatial resolution. The coefficients on Steller sea lion interactions are all positive and statistically significant except for MNLH, where it is negative and statistically significant. Given that the Steller sea lion interactions tend to occur in high-valued sites, this indicates that only the MNLH model is able to decouple the effect on latent site-specific heterogeneity, captured by the ASCs. This is consistent with the findings in Hicks and Schnier (2010). As will become more evident in the upcoming section, the effect of the Steller sea lion interactions is further reduced when the spatial resolution increases.
Focusing on the Cholesky factors indicates that all the sites have a lower variance than the second site, which is contained in Region 5410. This is consistent with the results in the six-site model as well. However, both Site 1 (the one normalized to zero) and Site 3 are also contained in 5410, which indicates that the Site 2 region (the central region of 5410) is more heterogeneous than the other two regions. Comparing the implied variance-covariance matrix for the random utility model using SMNP and SMNPH, 
, indicates that they are similar but not as similar as within the six-site model. This indicates that using the spatial autoregressive parameter to induce the spatial heterogeneity and autocorrelation within SMNP is not as effective as it was within the six-site model. Lastly, the results indicate that spatial autocorrelation is present, as p is statistically significant in both the SMNP and SMNPH models.
Twenty-Four-Site Model Results
The 24-site model generates results that are very similar to those of the 12-site model. The distance and missing dummy variable coefficients are always negative and statistically significant, and the revenue coefficient is positive and statistically significant across all models. The coefficients on the Steller sea lion interactions are different than those observed in the 12-site model. Within the 24-site model these coefficients are either negative and statistically significant or not significant at all. This indicates that when effort is measured at this spatial resolution, fishermen are actually deterred from fishing in locations where a higher number of Steller sea lion interactions occur. The results from the 6-site and 12-site models indicate the opposite. This illustrates the impact that spatial aggregation has on the behavioral model, as aggregating locations reduces the microvariation in the data, generating different implied behavioral responses.
Evidence of enhanced microvariation in the data can be found looking at the estimated Cholesky factors within the MNPH and SMNPH models, as there is a high degree of heterogeneity across the 24 sites varying between approximately 2 and 0.5 in the MNPH and 1.6 and 0.5 in the SMNPH models. Therefore, there is a sizeable degree of variation in the site-specific unobservables across the 24 sites. Restricting this variation to be constant across all sites pushes the Cholesky factors toward approximately 1 in the MNP model and 0.9 in the SMNP model, substantially reducing the implied variation in the unobservables. This reduced variation is also apparent when we compare the implied variance-covariance matrix for SMNP and SMNPH, 
, as the cross-correlation elements in SMNP are substantially smaller than the SMNPH results. Presumably these combined factors will compromise the performance of these estimators relative to the more flexible specifications of MNPH and SMNPH. This will be investigated in the following section. This said, the spatial correlation coefficient, ρ, is positive and statistically significant in both the SMNP and SMNPH models, indicating that spatial autocorrelation is present in the data.
Summary of Model Results
Summarizing the results from the three different spatial resolutions (18 models total) indicates that the distance traveled between sites is one the strongest predictors of spatial choice, as it is statistically significant and negative across all the models. Secondly, the expected revenues within a site predominantly have a positive effect on site selection. However, in the six-site model they are statistically insignificant and in some cases negative and statistically significant. This is a result of two factors: (1) the removal of the first implied haul from each constructed trip and (2) the spatial aggregation of effort into large regions that washes out some of the microbehavior observed within the data set.18 Removing the second element by expanding the spatial resolution reverses this phenomenon.
The coefficients on Steller sea lion interactions and the missing dummy variable were also impacted by the spatial aggregation of the data. At the smaller spatial resolution the Steller sea lion coefficients were predominantly positive and statistically significant, whereas at the larger spatial resolution they were either statistically insignificant or significant and negative. Similarly, at the smallest spatial resolution the missing dummy variable was positive and statistically significant, whereas at the smaller spatial resolution it flipped sign and was predominantly negative and statistically significant. These results should generate some caution for applied researchers because they illustrate that the parameter estimates obtained can be a construct of the data aggregation assumptions latent in the model. As was the case in our model, reducing the spatial resolution of the data also reduced the microvariation in the data that generated parameter estimates that did not accord with our a priori expectations. However, increasing the spatial resolution increased the observed microvariation, facilitating a better empirical model that accorded with our a priori expectations.
Spatial heterogeneity and autocorrelation are also prevalent within our data set. Although not illustrated within the empirical results, the ASCs estimated using the MNLH model were very heterogeneous, implying a high degree of spatial heterogeneity in the unobservables. All of the MNPH and SMNPH models estimated generated statistically significant differences in the site-specific Cholesky factors, illustrating that the variance in the site-specific unobservables is spatially heterogeneous. The spatial autocorrelation coefficients, ρ’s, in all the SMNP and SMNPH models were statistically significant, indicating that autocorrelation is present regardless of the spatial resolution of the data set. This said, comparing the implied variance-covariance matrix between the SMNP and SMNPH models indicated that increasing the spatial resolution of the data reduced the ability of the SMNP model to approximate the implied variance-covariance matrix generated from the SMNPH models. This result is primarily driven by the fact that the variance in the unobservables becomes more spatially heterogeneous as the spatial resolution gets smaller, thereby reducing the ability of the sole variance parameter combined with the spatial weighting matrix to approximate the full distribution of spatial unobservables obtained in the SMNPH model.
Model Comparisons
The discrete choice models estimated generate a probability surface that assigns a probability mass to each spatial alternative. Given that the likelihood functions used in the MNL and MNP (as well as SMNP) models are based on different error distributions, it is not possible to compare the log-likelihood values across the MNL and MNP (SMNP) models to determine which one provides a better statistical fit for the data. It is possible to use the log-likelihood to test across the different MNP and SMNP models, but given that one of the purposes of this research is to test these models against the MNL model, alternative performance metrics will be used. The first metric used is the model’s within-sample predictive accuracy. The second metric is each model’s out-of-sample predictive accuracy, and the third metric is the implied compensating variation resulting from a hypothetical site closure.
For the within-sample and out-of-sample predictive accuracies we utilize a measure of the model’s bias and mean squared error for comparison. The bias for each model is calculated as
[16]where Ii is a binary indicator for the observed choice within the data set, 
is the estimated probability for the spatial choice observed, and N is the total number of observations. The naïve (equal probability) predictions generate a bias measure of 0.833 for the 6-site model, 0.917 for the 12-site model, and 0.958 for the 24-site model.19 The mean squared error is calculated as
[17]where all variables are as defined earlier. The naïve predictions generate a mean squared error of 0.694 for the 6-site model, 0.840 for the 12-site model, and 0.918 for the 24-site model. The final performance metric, predictive accuracy, is calculated by selecting the choice alternative with the largest predicted probability and determining whether or not it was the actual spatial alternative selected in the data set. This metric is calculated as the percentage of correct predictions observed within the data set.
All three of these metrics are calculated for the years 2002–2005, within-sample predictions, and for the year 2006, out-of-sample predictions. The closer the bias and mean squared error calculations are to zero, as well as the greater the percentage of correct predictions, the better the model performs. To generate 95% confidence intervals for each metric, the Krinsky and Robb (1986, 1990) method was utilized with 120 draws from the parameter distribution.20
Within Sample Predictive Accuracy
The within-sample predictions for the 6- site, 12-site, and 24-site models are illustrated in Table 6. The models with the lowest bias, mean squared error, and predictive accuracies are bolded to differentiate them from the other results. For all the models estimated, the MNLH model possesses the lowest bias and mean squared error, with the MNL model following as a close second in many cases. In fact, there are only a few situations when the bias and mean squared error of the other four other models (MNP, MNPH, SMNP, and SMNPH) approach the bias and mean squared error of the MNL and MNLH models. This indicates that the probability surface estimated by these later models is flatter (i.e., assigns a higher probability mass to the nonvisited sites) than the MNL and MNLH models.
Within-Sample Predictive Accuracy of the Models: Average Values with 95% Confidence Intervals in Parentheses
The predictive accuracies within the six-site model illustrate that the MNLH model far outperforms the other models. The SMNPH and SMNP models follow this model, and the conventional MNL model actually outperforms the MNP and MNPH models. Comparing the models that do not incorporate spatial heterogeneity with the ones that do illustrates that accounting for spatial heterogeneity, either through the use of ASCs or site-specific Cholesky factors, increases the predictive accuracy of the model and also decreases the bias and the mean squared error. This relative increase is greater for the models based on the MNL framework than those based on the MNP framework. In the MNP framework the inclusion of spatial autocorrelation actually provides the largest increase in predictive accuracy relative to models that do not include spatial autocorrelation.
When our definition of space is reduced to the 12-site spatial resolution, the benefits of accounting for spatial heterogeneity and autocorrelation become a bit muted relative to the 6-site model. Furthermore, although the MNLH model does possess the lowest bias and mean squared error, it is actually marginally outperformed by all the other models when it comes to predictive accuracy. In fact, the predictive accuracy of the SMNPH model outperforms all the other models estimated. However, the confidence intervals across all the models, except the MNLH model, indicate that their relative performances are nearly identical and that the performance of the MNLH model is only marginally lower than the others. The lower predictive accuracy for the MNLH model is a result of the time-invariant nature of the ASCs. The contraction mapping used to estimate the ASCs is based on the aggregated site-specific probabilities over the entire data set. Therefore, although the ASCs do a very good job of predicting effort at this aggregate resolution, they do not perform as well at the microlevel.21 The fact that the bias and mean squared error are still lower than for the other models estimated, the predictive accuracy result indicates that when the MNLH model does not predict the correct discrete choice it misses the mark by only a small margin. Therefore, we cannot conclude that any one of the models at the 12-site spatial resolution is a clearly superior model.
At the 24-site spatial resolution, the benefits of accounting for spatial heterogeneity and autocorrelation become more evident. Although the MNL model does perform relatively well on average, the predictive accuracy of MNLH and SMNPH outperform the MNL model (results for MNPH model are nearly identical). However, the MNL does clearly outperform the MNP and SMNP models. The relative advantages of using the different models become more apparent when we compare the models that do not account for spatial heterogeneity with those that do account for spatial heterogeneity. Unilaterally accounting for spatial heterogeneity increases the predictive accuracy of the model and also reduces the bias and the mean squared error. In addition, when we compare the gains from accounting for spatial autocorrelation versus spatial heterogeneity we see that spatial heterogeneity provides the largest relative gains in model performance. This illustrates the benefits of using the MNLH, MNPH, or SMNPH models relative to their non-spatially heterogeneous counterparts.
Overall the within-sample predictions illustrate that accounting for spatial heterogeneity provides the largest relative gains in model performance. Furthermore, when the spatial resolution of the data is large (i.e., six sites), there exist some sizeable gains when accounting for spatial autocorrelation. This illustrates that the spatial independence assumption is not accurate and supports the conjecture that there exists a large degree of spatially omitted variables at this resolution of the data. At the 12-site model the gains from reducing this bias become more muted; however, the benefits of accounting for spatial heterogeneity are still present. Finally at the 24-site resolution the clear gains in model performance appear to be coming more so from accounting for spatial heterogeneity than spatial autocorrelation.
Out-of-Sample Predictive Accuracy
Although the within-sample predictive accuracies revealed that accounting for spatial heterogeneity and autocorrelation increases the statistical power of the models (granted, marginally in some cases), it is also important that these models generate reliable behavioral predictions, as these models may be used to guide policy. The out-of-sample metrics of performance are illustrated in Table 7, with the best performing models again highlighted in bold. With the exception of the fact that the bias and the mean squared error for the MNL model are actually lower than for all the other models at the 12-site resolution, the results are strikingly similar to the within-sample predictive accuracies. At the six-site spatial resolution the clear winner is again the MNLH model, followed by the SMNPH and SMNP models, clearly indicating that accounting for spatial heterogeneity increases the model performance. The 12-site model results are less pronounced, with all the models performing relatively well except for the MNLH model. This later finding presumably is a result of the ASCs that are capturing latent spatial information present in the 2002–2005 period that was not as present in 2006, hence the poorer predictions.
Out-of-Sample Predictive Accuracy of the Models: Average Values with 95% Confidence Intervals in Parentheses
In the 24-site model the out-of-sample prediction results are also very similar to the within-sample predictions. The MNLH model is the clear winner, and the MNPH and SMNPH models provide a sizeable gain in model performance relative to their non-spatially heterogeneous counterparts. In fact, the gains from accounting for spatial heterogeneity within the 24-site model are larger than those observed within sample. In addition, accounting for spatial autocorrelation actually reduces the out-of-sample performance. This indicates that the spatial correlation in unobservables is not consistent in our within-sample time period (2002–2005) and out-of-sample time period (2006). These results may be an artifact of the specific empirical application, but they do illustrate, yet again, the benefits of accounting for spatial heterogeneity.
Welfare Comparisons
To compare welfare estimates across the models Site 5431 was removed from a fisherman’s choice set.22 Given that the coefficient on the revenues in the 6-site model is predominately insignificant or does not accord with our a priori expectations, we focus on only the 12-site and 24-site model results. In both of these models all sites contained within 5431 were closed in order to facilitate comparison across the different spatial resolutions. Denoting X1 as the fully open choice set and X0 as the closed choice set, the compensating variation (CV) for the hypothetical site closure can be calculated as
[18]Estimating VMAX( · ) for MNL Models 1 and 2 is straightforward because
[19]where γ is Euler’s constant and Nj denotes the number of sites open in the fully open (j = 1) and closed choice sets (j = 0), respectively. For the MNP and SMNP models we utilize a frequency simulator (Chen and Cosslett 1998) to estimate VMAX( · ), which can be expressed as
[20]I[·] is an indicator function determining whether or not site i possesses the maximum site-specific utility measure conditional on the error vector draw 
We set R = 100 in our frequency simulator. The compensating variation estimates for the hypothetical closure are contained in Table 8.
Welfare Estimates Resulting from a Hypothetical Site Closure: 95% Confidence Intervals in Parentheses
The results from the 12-site model indicate that the MNL model generates the highest welfare estimates, whereas the 95% confidence intervals for the other five models are all overlapping, indicating that they are statistically similar. In fact, the welfare estimates for the MNL model are over twice those obtained from any of the other models. This is a striking result given that each of these models performs very similarly with respect to within-sample and out-of-sample performance measures. These results highlight the importance of not only accounting for spatial heterogeneity and autocorrelation. They also draw attention to fact that the MNL model may generate sizably different welfare estimates relative to a simple MNP model.
The 24-site model results indicate that the MNL and the MNP models generate very similar welfare estimates. The 95% confidence interval for the MNP completely covers that of the MNL model. Furthermore, our welfare estimates for the MNL and MNP models are all greater than the other four models estimated. This indicates that a failure to account for spatial heterogeneity and autocorrelation will increase our welfare estimates. On average, accounting for spatial heterogeneity in the MNL model, using the MNLH model, reduces the welfare estimate by 47%. In the MNP models, accounting for spatial heterogeneity alone reduces the welfare estimates by 12% (MNP versus MNPH) and accounting for both spatial heterogeneity and autocorrelation reduces it by 22% (MNP versus SMNPH). Therefore, the benefits from accounting for spatial heterogeneity are more sizeable in the MNL model. Furthermore, it is important to note that the confidence intervals for the MNP models are substantially larger than those observed in the MNL models. This indicates that the welfare estimates, although lower on average than the MNL model, are less precise in these models.
V. Conclusion
This research illustrates different empirical models that can be used to account for spatial heterogeneity and autocorrelation in discrete choice models. It also illustrates a methodology that may be used to account for both simultaneously, as well as estimate the off-diagonal variance-covariance parameters via parameterization of the spatial weighting matrix. Testing these models with a data set of spatial discrete choices using three alternative definitions of “space” illustrates that it is extremely important to account for spatial heterogeneity in the unobserved portion of spatial utility. In fact, the results indicate that using an MNL model with ASCs (MNLH) provides a very good control for spatial heterogeneity in discrete choice models. Furthermore, additional gains can be obtained when accounting for spatial autocorrelation in the MNP models. Therefore, should a researcher believe that the spatial independence assumption traditionally assumed in the MNL is restrictive, the MNPH and SMNPH models may provide a useful model for estimation.
However, there are some inherent limitations to using the MNPH and SMNPH models. They are not nearly as computationally tractable as the MNLH model. Because the probabilities must be simulated prior to maximizing the likelihood function using the GHK simulator, the computer time required to estimate the models increases substantially as the number of spatial sites increases.23 Furthermore, failure to account for spatial heterogeneity, captured by the MNP and SMNP models, generates results that are in some cases inferior to those of the MNL model. However, another alternative that may yield results that are superior to those of the MNL model would be to impose more structure on the error correlations in the MNPH and SMNPH models to reduce the number of parameters that are estimated. For instance, a spatial area could be divided into separate, larger spatial regions with each possessing identical error distributions combined with the spatial autocorrelation parameter. Undoubtedly the precise nature of the parameter restrictions, as well as the relative advantages of the different empirical models, will vary by the empirical application, but hopefully this research will motivate other researchers to explore the utilization of these methods and to ensure that they are taking sufficient steps to control for spatial heterogeneity and autocorrelation, where applicable, in their empirical estimates.
Footnotes
The authors are, respectively, associate professor, Department of Economics, Andrew Young School of Policy Studies, Georgia State University, Atlanta; and economist and program manager, Alaska Fisheries Science Center, U.S. National Marine Fisheries Service, Seattle, Washington.
↵1 Bockstael and Opaluch (1983) are credited with the first use of a discrete choice model in the fisheries literature in their study of effort supply decisions in the New England groundfish fishery.
↵2 The time subscripts for the site-specific utilities have been removed to simplify the model. For the purposes of this analysis, it is assumed that there is no serial correlation in the unobservables.
↵3 Other methods can be used to normalize the model of scale and level; see discussion by Train (2003). This transformation implicitly assumes that
and , which is necessary to facilitate estimation (Train 2003).↵4 A more detailed discussion of the simulated maximum likelihood routine is provided by Bolduc (1999) and Train (2003). Train (2003) describes in detail the steps required to utilize the GHK simulator.
↵5 Lower levels of K cut the estimation time substantially. Setting it equal to 50 roughly halved the estimation time, but levels lower than 50 generated unstable parameter estimates when using different random draws from the uniform distribution. In our analysis when K exceeded 50, the parameter estimates became more stable for different sets of random draws.
↵6 Hicks and Schnier (2010) investigate the endogenous formation of consideration sets invoked by this spatial regulation. We will not endeavor to add this dimension of complexity to the research to simplify the analysis. However, only those sites available to a vessel (as dictated by the spatial lottery) at a given point in time are used in the analysis.
↵7 Many of the sites included in this analysis are smaller than this spatial resolution because they are clumped around the islands within the BSAI and further subdivided for management purposes.
↵8 Vessel behavior residing outside of NMFS management zones 541, 542, and 543 were assigned a site value of 5000. Therefore, the six sites used in the larger spatial resolution are 5000, 5410, 5420, 5421, 5430, and 5431.
↵9 There are actually 49 distinct ADF&G zones contained in the data set, but many of them are visited only a few times in the analysis. Therefore, to simplify the analysis, sites that were visited less than 0.5% of the total hauls were aggregated into the nearest spatial neighbor.
↵10 Incorporating the first haul in the data set generates qualitatively similar results to those illustrated in the paper. These results can be obtained from the authors.
↵11 To control for exceptionally short and long trips, any trip that was less than 10 hauls or greater than 120 was not included in the data set.
↵12 As mentioned earlier a 30-day “burn-in” period was utilized within the data set. Lagged data from 2001 was not used because the regional TACs changed from 2001 to 2002, thus biasing the expected revenue calculations. Furthermore, all revenues were inflated to 2006 dollars assuming 3% annual interest.
↵13 The mean reported latitude and longitude were used because fishermen often cluster near shore, and using the centroids of the locations would bias the distance parameter.
↵14 More specifically, in our preliminary investigations we experimented with using an inverse distance, inverse distance squared (the model selected), and an inverse distance cubed spatial weighting matrix.
↵15 Although the log-likelihood value was the greatest using distance squared, our preliminary analysis indicated that the qualitative results reported within the paper do not depend on the selection of the distance weight. Furthermore, inverse distance squared is a common metric utilized in the spatial econometric literature.
↵16 To investigate the sensitivity of our results to this assumption we also experimented with using spatial weighting matrices that had the rows and columns for closed sites removed and rescaled. The results, as well as the model comparisons, from these regressions were very similar to those reported in the paper.
↵17 To investigate the robustness of this argument we also estimated the model including the first haul of the trip, assuming that the vessel originates in Dutch Harbor, Alaska (the busiest port within the area). For all of the models estimated, the coefficient on expected revenues was positive and statistically significant and the coefficient on the missing dummy variable was predominately insignificant. We elected to not use these results because of the potential measurement error they induce; however, these results are obtainable form the authors upon request.
↵18 Not removing the first haul of our constructed trips causes the revenue coefficients to become positive and statistically significant. Although we removed the first haul because we do not precisely know the true first haul of the trip, thus generating bias, there does have to be a first haul. Therefore, we can assume that these generalizations will hold if we know the “true” first haul of the trip.
↵19 The naïve expectations assume an equal probability of site visitation for all sites. Therefore, in the 6-site model the predicted bias is (1 - 1/6) = 5/6, for the 12-site model it is (1 - 1/12) = 11/12, and in the 24-site model it is (1 - 1/24) = 23/24. This information is provided as a baseline for comparison with the other model results.
↵20 Although 120 draws is smaller than the conventional number of draws utilized in the literature, it was selected due to the computational requirements of the multinomial probit models. Running 120 draws took approximately 10 hours for each MNP simulation.
↵21 Using time-varying ASCs would alleviate the problem; however, this is beyond the scope of this research effort. See Hicks, Horrace, and Schnier (2009) for an application that utilizes time varying ASCs.
↵22 It is important to note that our closure of site 5431 is purely hypothetical and in no way reflects the current management bodies’ interest in closing this region.
↵23 The 6-site SMNP model took approximately 15 minutes to estimate on a dual-quad core Mac machine, whereas the 24-site model took approximately 1 day to estimate.
