IntroductionTop

The origin of hedonic regression lies in valuing land of agricultural use (Haas, 1922) which, at the end of the 20th century and the start of the present century and with computers, has been well applied to value land worldwide (Xu et al., 1993; Shi et al., 1997; Maddison, 2000), and evidently in Spain (Caballer, 1973; Calatrava & Cañero, 2000; García & Grande, 2003; Gracia et al., 2004; Caballer & Guadalajara, 2005). In all these works, valuing models has been estimated by Ordinary Least Squared (OLS). However, spatial data, e.g. land values, present two properties that make meeting requirements and fulfilling the hypothesis of hedonic regression estimated by OLS difficult (Guadalajara, 2018): (1) geographic entities are spatially autocorrelated, (2) models are not stationary and perform differently in distinct study areas.

As a result, the OLS multiple regression estimations for the βi coefficients of the independent variables will be most probably biased and inconsistent, and will also invalidate standard regression diagnostic tests through misstated standard errors (Kim et al., 2003).

Autocorrelation, association or spatial dependence refers to the concentration or dispersion of the values of a variable (land prices in our case) in a land or geographic space. This implies that the value of a variable is conditioned by the value that this same variable takes in one neighboring region or in several. According to the first law of geography of Tobler (1970), “Everything is related to everything else, but near things are more related than distant things”.

The first spatial studies were done after hedonic models at the end of the 20th century (Can, 1992; Pace & Gilley, 1997; Dubin, 1998) thanks not only to geographic information being implemented and access to big databases gained, but also to Geographic Information Systems (GIS) and software being de­veloped to analyze spatial data. In these GIS, data are geo-referenced by latitude and longitude, or by Universal Transverse Mercator (UTM) X Y coordinates (Guadalajara, 2018). Spatial regression models applied to land valuing have been well developed in the present century. Generally speaking, the most widely used spatial models are the spatial lag model (SLM) and the spatial error model (SEM), and both are applied to correct spatial autocorrelation. SLM includes a spatially lagged dependent variable, while SEM includes the spatial dependence of the error term. Some examples of these spatial models that have been applied to land valuing are those by: Patton & McErlean (2003) in Northern Ireland; Huang et al. (2006) in the USA; Seo (2008) in South America; Maddison (2009) in the UK; Mallios et al. (2009) in Greece; Zygmunt & Gluszak (2015) in Poland; Uberti et al. (2018) in Brazil. They all indicate the need to consider GIS because hedonic models do not fulfil the OLS hypothesis.

The models obtained in the above-cited works include two main categories of explanatory variables: internal and external in relation to property. We cite the following internal variables: (1) irrigation availability: this variable is considered a dummy variable and takes a positive sign in relation to not only the land unit price logarithm found in the work by Mallios et al. (2009), but also to the land unit price in the work by Demetriou (2016), insofar as irrigated land increases the land unit value; (2) plot size: in some cases this variable is considered in its original form (Patton & McErlean, 2003; Maddison, 2009; Zygmunt & Gluszak, 2015; Demetriou, 2016) and in a logarithmic form in others (Huang et al., 2006; Mallios et al., 2009), but it always takes a negative sign in relation to the unit price logarithm (Huang et al., 2006; Maddison, 2009; Mallios et al., 2009) or the land unit price (Patton & McErlean, 2003; Zygmunt & Gluszak, 2015; Demetriou, 2016). This indicates that land unit values lower with plot size; that is, the total price of plots does not linearly in­crea­se with surface; (3) topography: Demetriou (2016) obtained a negative relation between the land unit price and plot slope as so: the steeper the slope, the lower the unit price; (4) altitude: Mallios et al. (2009) obtained a positive relation between unit price and land rise, and both logarithmically insofar as lands at higher altitude are more expensive; (5) the arable land classification influences the land unit price, which takes a negative sign in some works (Patton & McErlean, 2003; Maddison, 2009) and a positive one in another study (Huang et al., 2006), depending on how the arable land classification is quantified.

We cite the following external variables, among others, and numerous variables controlling for locational effects: (1) distance to reference places like: (a) distance from residential zones, which always has a negative effect on the land price (Bastian et al., 2002; Patton & McErlean, 2003; Huang et al., 2006; Maddison, 2009; Mallios et al., 2009; Zygmunt & Gluszak, 2015); (b) presence of sea: both Demetriou (2016) and Mallios et al. (2009) consider sea views and the distance from the sea logarithm, respectively, with a positive sign for the land unit price, which implies an opposite effect in each case; (c) distance to the nearest main road appears in the models of Mallios et al. (2009) with a negative sign for the land unit price logarith. However, Uberti et al. (2018) and Demetriou (2016) report a positive relation between access to plots and the unit price, which means the same in all three cases: better accessibility to plots increases their price; (d) distance to forest negatively impacts land unit prices (Zygmunt & Gluszak, 2015); (2) population density and personal income per capita (Huang et al., 2006) influence prices insofar as the land values logarithm increases with population density and personal income per capita, and both logarithmically; (3) environmental amenities, like fishing (Bastian et al., 2002), positively influence the land unit price logarithm.

In Spain, GIS have been used to model the location factor set out in the Spanish Land Act (Marqués-Pérez et al., 2018). Although the spatial correlation of land values has been demonstrated (Segura & Marqués, 2018), no spatial models have been obtained to explain arable land values, only for house values (Militino et al., 2004; Taltavull et al., 2016; Guadalajara & López, 2018).

Consequently, the objective of the present work was to obtain spatial models to value land used for agriculture by distinguishing among the uses that collect the spatial autocorrelation of land values, and to improve the results obtained with conventional models. At the same time, the intention of using these models was to determine and quantify the factors explaining land prices. The data employed to obtain these models were the mean prices per municipality and per land use type in the Spanish Autonomous Community of Aragón (SACA).

Material and methodsTop

Data

The employed information source was the website of the Aragonese Statistics Institute (IAEST, in Spanish) of the SACA. The SACA covers 47,720 km2 and is divided into three provinces: Huesca to the north, Zaragoza in the center and Teruel to the south. The following information was collected for the 741 municipalities in the SACA in June 2018: (1) mean land price (€/ha) per land use in 2017; (2) the internal characteristic in relation to property: geographic coordinates (longitude and latitude, time zone of the UTM projection and X UTM and Y UTM); agricultural characteristics (usable agricultural area (UAA [ha]), irrigatable area in relation to UAA in percentages and number of plots on rustic land); orographic characteristics (altitude [m]); (3) external characteristics in relation to property: demographic characteristics (population size; population's mean age; birth rate; death rate); services (number of compulsory secondary education centres (CSEC)); economic characteristics (cadastral value of rustic land in thousands of euros of the whole municipality and gross per capita income in 2014 (euros per person and year) in seven intervals: < 6000; 6000–7999; 8000–9999; 10000–11999; 12000-15999; 16000-17999; ≥18000). With the above information, the mean plot size of each municipality (UAA/no. plots) and population density (population/UAA) were calculated to include them in the study, which falls in line with previous works.

The land use types in the IAEST were: almond trees (non-irrigated and irrigated), arable land (non-irrigated and irrigated), olive groves (non-irrigated and irrigated), vineyards (non-irrigated and irrigated), meadows (non-irrigated and irrigated), irrigated fruit trees, orchards, wasteland, pinewoods and riverside trees. The surface of each land use type was calculated using the 2018 Surface Areas and Crop Yields Survey, with the results summarized by SACA (Spanish Ministry of Agriculture, Fishing and Food: www.mapa.gob.es).

Similarly to other works that have considered distance to places of interest, the location of municipalities in some of the 18 nature reserves in the SACA was also contemplated, which are listed on this website: www.aragon.es/-/red-de-espacios-naturales-protegidos. In all, 6686 observations made up the analyzed sample.

Maps were created with mean prices and for each land use per municipality, shown in quantile inter­vals using ArcGIS Pro 2.2.0 (©2018 Esri Inc.). The UTM projection system and the reference ETRS89 geodesic system were used, time zone 30N, in which the municipalities forming part of time zone 31 were included. The shapefile employed in the areas corresponding to municipality limits (recintos_muni­cipales_inspire_peninbal_etrs89.shp, type: ‘Poly­gon’, uncertainty range = 40m, download date 23 July 2018) was obtained from the Centro de Descargas del Centro Nacional de Información Geográfica del Instituto Geográfico Nacional (the Downloading Center of the National Geographic Information Center of the Spanish National Geographic Institute; the Spanish Ministry of Development, the Spanish Government, www.ign.es). This center lists the Aragonese municipalities that form part of time zones 30 and 31.

Testing for spatial effects

In order to test the presence of spatial effects in the data about mean prices per municipality, spatial weights wij between municipalities were calculated. Weights represent the geographical relationship bet­ween locations i and j. Several methods are available that construct spatial weights: contiguity (Huang et al., 2006), k-Nearest Neighbor (Zygmunt & Gluszak, 2015; Uberti et al., 2018) and distance (Patton & McErlean, 2003; Maddison, 2009; Zygmunt & Gluszak, 2015; Uberti et al., 2018). As spatial information comes as geographical coordinates (point data), this work intended to build weights by considering the distance among municipalities, as most authors have done, by using the values X UTM and Y UTM. Weights were calculated in two ways: by taking the inverse of Euclidean distance squared, as Patton & McErlean (2003) did, and with the inverse of Euclidean distance, as Maddison (2009) did, to select that which would provide the most compelling evidence for spatial dependence. For all land uses, a minimum threshold distance was considered so that all the municipalities had at least one neighbor which, at the same time, would be the maximum permitted distance to consi­der a municipality a neighbor. Spatial weights matrix W = [wij] contains weights between each pair of all observations (municipalities) and is a non-negative m×m matrix. Matrix elements cannot be their own neighbors insofar as the matrix's diagonal line is composed of zeros. The weight matrix was standardized in such a way that the sum of the weights in each row equaled 1.

Having defined the spatial weights and the spatial weights matrix, the global Moran’s Index (I) test statistics (Moran, 1950) was used, which is the most popular statistics to measure spatial association, whose value varies between -1 (perfect dispersion) and 1 (perfect correlation). A value of 0 indicates a null correlation or a random spatial pattern, and the nearer it comes to 1, the higher the spatial correlation.

Regression models

The methodology used to obtain land valuing models was hedonic regression models. First an estimation by OLS was done. The basic linear hedonic model, using the log-linear model (Pace & Gilley, 1997; Bastian et al., 2002; Maddison, 2009; Mallios et al., 2009; Zygmunt & Gluszak, 2015), is given by:

where dependent variable Yi is an m×1 vector of the mean land value for each municipality (m is the number of municipalities); α is the constant term; Xij is a m×n matrix of the independent variables (n is the number of explanatory variables); βj is a n vector of the estimated coefficients of the independent variables, and εi is an m×1 error term. Independent variables can be quantitative or dummy, and quantitative variables can come in their original form or be transformed into a logarithm. If variables are quantitative, each coefficient βj represents the elasticity of demand for this speci­fic characteristic (Gujarati, 2003). If characteristic Xij comes in its original form, when Xij varies by 1 unit, then Y varies by βj·100%. on average. If the characteristic comes in a logarithmic form, when Xij varies by 1%, then Y varies by βj%. on average. If a characteristic comes as a dummy variable, coefficient βj provides the exact percentage difference between the inclusion of the characteristic or not as expβj -1 (Mallios et al., 2009).

Eleven models were obtained, one for the set of lands in the SACA and 10 other models for all ten conside­red land uses. Initially, all the aforementioned variables from the IAEST were included as independent variables. The model for the set of lands in the SACA also included nine dummy variables relating to land use, which took a value of 1 if they were related to the land use in question, and 0 otherwise. To distinguish between non-irrigated and irrigated land uses, another dummy variable was included, namely “Irrigation”, which took a value of 1 if it was an irrigated crop, and 0 otherwise. The dummy variable “Nature reserves” was also contemplated, which took a value of 1 if the municipality was located in a nature reserve, and 0 otherwise. A municipality's altitude was considered in km and the number of plots in thousands.

Quantitative variables: cadastral value and population were considered in two ways: in their original form and in their transformed logarithmic form. The municipality's income took values from 1 to 7, with 1 corresponding to the lowest income interval and 7 to the highest.

In order to begin the spatial regression analysis, the spatial autocorrelation in the OLS residuals was evaluated by Moran’s I test, done with the residuals to ensure that they were spatially random. The spatial matrix captures the spatial autocorrelation present in the residuals of the hedonic regression by OLS. The spatial regression model specification tests were obtained ba­sed on the Lagrange Multiplier (LM) of the dependent variable, LM-lag, and of the error, LM-error, and also in their robust versions. These tests allowed the problem of the specification of the spatial regression models to be solved. Thus we considered two spatial regression models to incorporate the spatial components into the OLS (Anselin, 1988):

The Spatial Lag Model (SLM) or the Spatial Autoregressive Model (SAR):

According to this model, a land value is considered to be autocorrelated in space. This model is formally written as:

where WlogYi is the spatially lagged dependent variable (additional regresor) and ρ is the spatial lag coefficient. WlogYi is defined as:

The spatially lagged dependent variable is interpreted as a weighted average of the neighboring land values.

The Spatial Error Model (SEM):

This model handles spatial dependence through the error term, and takes the following form:

where λ is the coefficient of the spatially correlated errors and Wui is the spatially lagged error term.

According to Anselin (1988), the estimation of models SLM and SEM cannot be done by OLS, but by using Maximum Likelihood (ML), which is based on the normality and independence hypotheses of the error term.

To evaluate the goodness-of-fit criteria, the deter­mination coefficient (R2), the log likelihood, Schwarz's criterion (SC) and Akaike’s information criterion (AIC) were used to test several functional forms for the hedonic price equation and the selected variables, and also the SLM and SEM models estimated by ML. SC, AIC and log likelihood are an appropriate measure for comparing non-nested models. Models with smaller AIC and SC are considered superior (Chi & Zhu, 2008). Conversely, the higher the log likelihood value, the better the model fits. The Likelihood Ratio Test (LRT) is a statistical test of the goodness-of-fit between two spatial models.

The procedure followed to select the variables was the stepwise method. A Student's t-test was done of the coefficients of all the independent variables to define the variables that were significant and had influenced the value. Those variables with non-significant coefficients were not included in the final model.

For the regression diagnostics, the collinearity or combination of the explanatory variables was determined by the condition index (CI), and was also analyzed by the variance inflation factor (VIF) of each explanatory variable. Gujarati (2003) indicates serious multicollinearity problems likely exist with condition index scores over 30 and recommends a lower VIF than 10 (rule of thumb threshold).

For the regression diagnostics, the Koenker-Bassett (K-B) and Jarque-Bera (J-B) statistics were used in the OLS model. If the K-B and J-B statistics are statistically significant, they indicate that the distribution of the residuals is not normal, and the OLS results will have an incorrect model specification (a key variable that is lacking in the model). The Breusch-Pagan (B-P) statistics was used to test all the regression models. If the B-P statistics is significant, it indicates that the model is not consistent. That is, the relations being modeled change in the study area (non-stationarity) or vary in relation to the magnitude of the variable that is to be foreseen (heteroscedasticity). The GeoDa software was used to obtain Moran’s I test, as were the OLS, SLM and SEM models with their statistics.

ResultsTop

Spatial effects

The number of municipalities for which information existed about prices for land uses, the mean, minimum and maximum price values, Moran’s I test corresponding to these prices, and the surface of each land use type are found in Table 1. The analyzed uses represented 65.35% of the SACA’s surface area, where non-irrigated arable land (25.47%) predominated, followed by pinewoods (16.40%). The following were not included because their price information was not available: non-irrigated fruit trees, scrubland, thickets and conifers, among others. This table also includes the threshold distance, that is considered to calculate the spatial weights, for which all the municipalities have at least one neighbor. The spatial weights were calculated with the inverse of Euclidean distance because it provided identical Moran’s I test values to the inverse squared.

As Table 1 shows, the mean price per municipality in the SACA ranged from a minimum of 120 €/ha for wasteland to a maximum of 33,640 €/ha for irrigated orchards, and the mean value was 4317 €/ha. High Moran’s I test values indicated that a high spatial correlation exists in the land prices for all land uses, except for irrigated meadows and irrigated lands with fruit trees, for which Moran’s I was only 0.0072 and 0.0783, respectively. The highest spatial correlation of land prices was obtained for non-irrigated meadows (0.9143), followed by riverside trees (0.8754), wasteland (0.8742) and non-irrigated arable land (0.8637).

The maps showing the mean values per municipality for each land use, represented in price intervals by quantiles, are found in Figures 1-4. They all confirmed that a high spatial correlation existed for land values, except for irrigated meadows and irrigated lands with fruit trees. For all land uses, the highest prices for non-irrigated land were obtained for the province of Huesca, for the irrigated lands in the Ebro Valley and to the east of Huesca. Conversely, the province of Teruel obtained the lowest prices and also the fewest different land uses.

Figure 1. The land value (€/ha) in the municipalities of the Spanish Autonomous Community of Aragón of (a) non-irrigated arable land, (b) irrigated arable land, (c) non-irrigated almond trees, and (d) irrigated almond trees.

Figure 2. The land value (€/ha) in the municipalities of the Spanish Autonomous Community of Aragón of (a) non-irrigated olive groves, (b) irrigated olive groves, (c) non-irrigated vineyards, and (d) irrigated vineyards.

Figure 3. The land value (€/ha) in the municipalities of the Spanish Autonomous Community of Aragón of (a) non-irrigated meadows, (b) irrigated meadows, (c) irrigated fruit trees, and (d) orchards.

Figure 4. The land value (€/ha) in the municipalities of the Spanish Autonomous Community of Aragón of (a) riverside trees (b) pinewoods, and (c) wasteland.

Regression models

Table 2 includes the OLS, SLM and SEM models for all the land uses in the studied SACA, where wasteland use is considered as control. Table 3 shows the OLS models that corresponded to each land use by grouping non-irrigated and irrigated in those land uses where both these possibilities were given.

Tables 2 and 3 show that the LM-lag and LM-error statistics were both significant and, therefore, the ro­bust versions of the statistics were taken into account. Both the robust and non-robust versions of the test statistics were significant, except for the robust LM-lag for pinewoods. Therefore, both spatial models were suitable for modeling land values for their different uses, including fruit trees and irrigated meadows. Nevertheless, following Anselin & Rey (1992), the results for LM-lag and LM-error shown in Tables 2 and 3 could indicate that the SEM was the most appropriate model to describe the land value of pinewoods, as well as meadows, irrigated lands with fruit trees, irrigated riverside trees and wastelands because the LM-error values were higher than the LM-lag values. Conversely, the SLM would be more appropriate for arable land, almond trees, olive groves, vineyards and orchards, and to also describe the set of all land uses.

The highest CI scores were 31.62 for lands with almond trees, followed by 31.48 for wastelands. The CI scores were always below 30 for all other land uses. As all the VIFs were below 3, all these diagnostics indicated that no multicollinearity existed in these models.

The normality of the residuals was not met, as the J-B test results revealed. So the null hypothesis of a normal error was rejected. The exceptions were wastelands and meadows, for which the J-B test was not significant, but was able to confirm the normality of the residuals.

Tables 4 and 5 respectively show model SLM and model SEM, which correspond to each land use. In order to select the best model, and in accordance with R2, the log likelihood, the AIC and SC, the spatial models were always superior to those obtained by traditional OLS. Considering the spatial models im­proved the goodness-of-fit measures. R2 was always higher in the spatial model, unlike OLS, especially in riverside trees (0.90 vs. 0.36), wasteland (0.86 vs. 0.38) and pinewoods (0.76 vs. 0.31). The same was true of the log likelihood, which increased in the spatial models, especially in the SEM models for the set of land uses (from -1252.11 to -204.18), riverside trees (from -3.55 to 598.16) and wasteland (from -144.55 to 327.88).

AIC and SC lowered in all the spatial models. AIC went from 2538.21 to 440.36 for the set of land uses, from 21.10 to -1188.33 for riverside trees and from 303.09 to -647.77 for wasteland. Respectively for the same uses, SC went from 2653.85 to 549.20, from 53.16 to -1175.32 and from 335.16 to -646.94.

The fact that all the spatial autoregressive terms were significant (ρ for SLM and λ for SEM), indicated substantial spatial effects across the municipalities in the SACA (Lee et al., 2016). According to Huang et al. (2006), spatial autoregressive estimate ρ, which ranged between 0.2567 for the model for the value of the set of land uses and 0.8962 for the value of the riverside trees, indicated that a 1% increase in the average land prices in nearby municipalities would increase the land prices in the observed municipality by 0.2567% and 0.8962%, respectively. The high positive ρ value indicated that the land value was strongly influenced by the values of neighboring lands. The ρ values were higher in the models for uses than in the model for set of land, which meant that this neighborhood effect was more marked on the land values for uses than on the land value in general. The same occurred with coefficient λ in relation to the correlation of the residuals, which was higher in the models for uses. This gave way to most of the coefficients of the other explanatory variables being higher in OLS than in spatial models because spatial coefficients ρ and λ collected part of the land values owing to the neighborhood effect.

The LRT result was significant in all the models obtained for the different uses, including the model for lands with pinewoods for which, as we have found before, the robust LM-lag was not significant and obtained similar values in both spatial models. Therefore, it was corroborated that the two spatial models were suitable for modeling land prices.

The homoscedasticity of the residuals was verified by the B-P test and the results indicated that the residuals had non-constant variance and heteroscedasticity problems in all the models, including the spatial models. Once again, the exception went to the model for meadows, where the homoscedasticity of the residuals was confirmed in the SLM model, and in the lands with olive groves.

The results presented in Tables 2-5 confirmed the importance of both agricultural and non-agricultural factors for determining land prices.

According to the coefficient values of the uses obtained for the model with the set of land uses (Table 2), the land value in relation to the wastelands value oscillated from 1.57-fold (exp0.9426-1), or 157% for lands with riverside trees, to 11.95-fold (exp2.5612-1) or 1195% for orchards.

In the models obtained for the set of land and for each use, the explanatory variables differed, but all the coefficients were significant at the 5 percent level at least, and the signs of the coefficients corresponded to a priori expectations. Irrigated land always took a positive sign because water availability always increases land productivity. According to the coefficient values of irrigation, in the model for the set of land uses in the SACA, which ranged from 1.1843 in the OLS model to 1.1609 in the SEM model, the difference between the price of irrigated land and that of non-irrigated land was 2.27-fold (exp1.1843-1) or 227% and 2.19-fold (exp1.1609-1) or 219%, depending on the model. For uses, the coefficient values ranged between 0.9055 for meadows in the SEM model and 2.3970 for vineyards in the OLS model. Hence the difference between the price of irrigated land and that of non-irrigated land ranged between 1.4732-fold (exp0.9055-1) or 147.32% for meadows and 9.9902-fold (exp2.3970-1) or 999.02% for vineyards.

Our results confirmed that land values increased with plot size and lowered with the number of plots in the municipality. Indeed for the land set, a 1 ha increase in plot size increased the land value by 6.49-7.88% (β = 0.0649 in SLM and β = 0.0788 in SEM). For land uses for fruit trees, a 1 ha increase in plot size in the OLS model increased the land price by 14.38% (β = 0.1438), while for riverside trees and according to SLM, a 1 ha increase in plot size increased the land price by only 1.09% (β = 0.0109).

In relation to the municipality's UAA, irrigated land areas only intervened in the model for vineyards and took a negative sign, while this characteristic did not appear in the model for the other models for each land use. Moreover, the UAA in its logarithmic form only appeared in the model for fruit trees and took a negative sign.

The municipality's income explained the mean price in the model for the land set and in the set of all land uses, except for vineyards and always with a positive sign. This was expected because it is indicative of a municipality's wealth, which tends to come with a higher land price. An increase in income within one interval gives way to a general increase in land of 4.77-5.63% (β = 0.0477 in SLM and β = 0.0563 in SEM). For land uses, this increase varied from 0.82% (β = 0.0082) for riverside trees according to SLM to 11.87% (β = 0.1187) for meadows according to OLS.

Another indication of a municipality's wealth is its cadastral value, which increases the land value for all land uses, except for meadows, fruit trees, orchards, pinewoods and riverside trees. A 1% increase in the cadastral value increased the land price by between 0.0251% (β = 0.0251) according to SLM and 0.0761% (β = 0.0761), according to OLS, and for arable land in both cases.

A bigger municipality population increased the land prices depending on the model for the set of land uses, and per use for arable land, irrigated fruit trees and orchards. Population density also increased the vineyard land value.

The municipality's mean age only influenced the vineyard land value and negatively so; i.e., the muni­cipalities with an older mean age obtained a lower vineyard land price. The death rate also had a negative influence, but only on the value in the set of land in the SACA.

A higher altitude lowered the land price for almond trees, irrigated fruit trees, orchards, wasteland and ri­verside trees, but the opposite occurred for lands with pinewoods and vineyards.

Finally, the location of a municipality in a nature reserve increased the land value in general, and for these uses in particular: almond trees, meadows, orchards, pinewoods, wasteland and riverside trees. According to the coefficient values of nature reserves, which ranged between 0.1046 and 0.1213 depending on the models for the set of land in the SACA, the difference between the price of land located in a nature reserve and that outside a nature reserved ranged from 0.1102-fold (exp0.1046-1) or 11.02% to 0.1289-fold (exp0.1213-1) or 12.89%. For uses, land values rose from 0.0344-fold (exp0.0339-1) or 3.44% for lands with riverside trees in SLM to 0.3531-fold (exp0.3024-1) or 35.31% for meadows in OLS for those municipalities located in a nature reserve.

DiscussionTop

This paper shows that spatial effects are significant on land values in the SACA, which coincides with previous studies conducted in other areas (Patton & McErlean, 2003; Huang et al., 2006; Seo, 2008; Maddison, 2009; Mallios et al., 2009; Zygmunt & Gluszak, 2015; Uberti et al., 2018).

The spatial correlation of land prices in Spain was confirmed, which falls in line with what Segura & Marqués (2018) obtained for mean prices per province. The Moran’s I test values herein obtained were higher for the analyzed uses: 0.8637 vs. 0.3475 for non-irrigated arable land, 0.7385 vs. 0.3755 for irrigated arable land and 0.9143 vs. 0.3075 for non-irrigated meadows. This could be due to the different reference levels employed, e.g., provincial vs municipal, as in our case.

Spatial models SLM and SEM proved better than OLS models for all the possible land uses, which also happened in the consulted studies. This indicates the need to develop spatial models to model land prices by implementing GIS. The LM-lag and LM-error statistics pointed out that SEM was slightly better than SLM in some models, which implies that the spatial effect was stronger on errors than on land prices. This could be due to some of the variables not being included in models, such as temperature, soil quality and precipitation. However, these data were not available for municipalities. This was corroborated by the significance of coefficient λ, which suggests that other explanatory variables may have been omitted from the models.

The R2 values obtained in the models developed herein were generally similar to those obtained by Huang et al. (2006), and were even higher than those reported in most of the consulted works: 0.60 in Bastian et al. (2002); 0.63 in Patton & McErlean (2003); 0.49 in Maddison (2009); 0.52 in Zygmunt & Gluszak (2015); 0.69 in Uberti et al. (2018).

Similarly to other works (Bastian et al., 2002; Patton & McErlean, 2003; Mallios et al., 2009; Demetriou, 2016; Guadalajara & López, 2018), it was not possible to eliminate the heteroscedasticity of the residuals in most of the models obtained for the land value in the SACA, as deduced from the B-P test results. Heteroscedasticity was eliminated only in olive groves and meadows, and lowered in all crops, except for vineyards and riverside trees, when spatial models were utilized. Apart from employing spatial models, a widely used resource to reduce heteroscedasticity is variables transformed into logarithms, which was done, but was not entirely successful. The inclusion of the municipality’s precipitation in the models could have lowered heteroscedasticity. Nonetheless, it is noteworthy that other consulted works (Huang et al., 2006; Seo, 2008; Maddison, 2009; Zygmunt & Gluszak, 2015; Uberti et al., 2018) did not indicate the result of either this test or the J-B test, which apparently suggests a problem in these models that needs to be solved. A literature review indicates that a joint remedy is lacking for these conditions when the nature of heteroscedasticity is unknown.

The multicollinearity condition number in the ob­tained models was lower than that indicated in other works, e.g.: 34.98 in Mallios et al. (2009) and 48.12 and 93.8 in Patton & McErlean (2003), which ratifies the importance and validity of the models developed in the present work.

The signs of the significant agricultural production variables met a priori expectations. Irrigation was al­ways positive, exactly as indicated by Bastian et al. (2002), Mallios et al. (2009) and Demetriou (2016). Irrigatable areas in relation to the municipality's UAA only influenced vineyards and took a negative sign. This could be due to a larger irrigatable surface area in relation to the total surface area, which could increase the supply of irrigated land and could lower its price. The municipality's mean age also influenced vineyards because land was demanded more in the municipalities with a younger population, which could have something more to do with the younger population's interest in producing wines.

Unlike other works (Huang et al., 2006; Maddison, 2009; Mallios et al., 2009; Zygmunt & Gluszak, 2015; Demetriou, 2016), land unit values increased with plot size. This ratio between unit values and plot size might depend on the characteristics of the crops in each country. In Spain, large surface areas mean mechanisation and lower crop costs. These lower land prices for smaller plot sizes are related with the findings reported by Levers et al. (2018) about which determining factors related to farm management, e.g. smaller field size, would contribute to describe agricultural abandonment patterns in Europe. The above authors' study indicated some areas in Spain, like Galicia and south Aragón, where the smaller the plot size, the more likely abandonment is.

Altitude also influenced values, which came over in the work by Mallios et al. (2009) but, in our case, it only had a positive effect on lands with pinewoods and vineyards. It negatively affected lands with almond trees, irrigated fruit trees, orchards and riverside trees, most certainly because these land uses are more sensitive to damage caused by low temperatures, which occurs more frequently at higher altitudes.

As maintained by Huang et al. (2006), land values increase with population density and personal per capita income. A denser population places more pressure on land use and leads to higher prices.

The influence of nature reserve locations on arable land prices was also shown and coincides with other works (Bastian et al., 2002) and also with the Spanish regulations (BOE, 2011).

The results of this study might be interesting for rural land management, the mass appraisal for the determining factor of market values, territorial taxation, and for actions to avoid land being aban­doned. One study limitation is the availability of the municipal data instead of data about plots, which did not allow us to include some specific plot characteristics, like plot shape (Zygmunt & Cluszak, 2015), plot slope (Demetriou, 2016), soil type, distance from the population center, etc., which can influence prices, as other studies have demonstrated. One future research line will consider the influence of proximity to communication routes (main roads, high-speed trains, etc.) and how they improve land prices, and to contemplate the protected designations of origin of some crops like wine.