# Fuzzy group fixed-effects estimation with spatial clustering

Roy Cequeti, Pierpaolo D'urso, Raffaele Mattera

<figure><img src="image_1.png" /></figure>

## To cite this version:

Roy Cequeti, Pierpaolo D'urso, Raffaele Mattera. Fuzzy group fixed-effects estimation with spatial clustering. AStA Advances in Statistical Analysis, 2025, 10.1007/s10182-024-00522-1. hal-05109570

HAL Id: hal-05109570

https://univ-angers.hal.science/hal-05109570v1

Submitted on 12 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Fuzzy group fixed-effects estimation with spatial clustering

Roy Cequeti¹,² · Pierpaolo D'Urso¹ · Raffaele Mattera¹,³

<figure><img src="image_4.png" /></figure>

Received: 17 March 2024 / Accepted: 4 December 2024

© Springer-Verlag GmbH Germany, part of Springer Nature 2025

## Abstract

The paper discusses the problem of estimating group heterogeneous fixed-effect panel data models under the assumption of fuzzy clustering, that is each unit belongs to all the clusters with a membership degree. To enhance spatial clustering, a spatio-temporal approach is considered. An iterative procedure, alternating panel data estimation and spatio-temporal clustering of the residuals, is proposed. The proposed method can be of relevance to researchers interested in using fuzzy group fixed-effect methods, but want to leverage spatial dimension for clustering units. Two empirical examples, the first on cigarette consumption in the US states and the second on non-life insurance demand in Italy, are presented to illustrate the performance of the proposed approach. The spatial fuzzy GFE model reveals important regional differences in both the US cigarette consumption and non-life insurance determinants in Italy. In the case of the US, we found a distinction in two main clusters, East and West. For the Italy provinces data, we find a distinction in North and South clusters. Regarding the regression results, for cigarette consumption data, different from the previous studies, we find that the smuggling effect is significant only in east regions, thus suggesting localised impacts of bootlegging. In the context of Italian non-life insurance demand, we find that while population density explains insurance consumption in northern provinces, the trust issues in the south explain the lower insurance demand.

**Keywords** Longitudinal data · Clusterwise regression · Grouped fixed effect · Fuzzy approach · Heterogeneous coefficient · Spatial econometrics

## 1 Introduction

A common starting point in regression analysis is the assumption that all statistical units have a uniform relationship between dependent and independent variables. However, real-world phenomena often exhibit complexities that challenge this assumption in practice. For example, there may be distinct groups of units with inherently different relationships. This problem, known in the statistical literature as

<footnote>Extended author information available on the last page of the article</footnote>
Simpson's paradox (Alin 2010), arises in empirical research when there are trends in subgroups of data that disappear or reverse when these groups are excluded.

Consider, for example, the problem of examining the factors influencing cigarette consumption in different regions. While aggregated results might point to common determinants (Baltagi and Levin 1992), a closer look at the data might reveal that these determinants have different effects on different demographic or geographic subgroups. Another illustrative example is the analysis of the determinants of the demand for insurance products (Gruber and Poterba 1994), where, for example, the factors influencing consumption among young adults may be very different from those affecting older populations, or where cultural factors common to groups of countries or regions may affect both the sign and the magnitude of the regression coefficients.

In standard cross-sectional regression, several clusterwise regression approaches have been proposed. Späth (1979) proposed one of the first clusterwise regression algorithms, alternating parameter estimation step and residuals based clustering, which was extended to fuzzy clustering setting by Hathaway and Bezdek (1993) and Coppi et al. (2006); D'Urso and Santoro (2006) and D'Urso et al. (2010), among the others. DeSarbo and Cron (1988) proposed a Gaussian mixture-based clustering regression approach, whose properties have been studied in more detail by Hennig (2000) and Brusco et al. (2008), while more recently, Sugasawa (2021) proposed a more general parametric approach for cross-sectional data. Di Mari et al. (2023) introduced a variable selection step in such a framework through LASSO-type shrinkage.

In the context of panel data regression, Simpson's paradox emphasises the importance of recognising and accounting for group-specific dynamics. Otherwise, conclusions may be misleading because aggregate results may mask significant variation within individual subpopulations. This issue has only recently been addressed in the panel data literature. To deal with it, advanced techniques such as group heterogeneous (also called clustered) panel data models (Bester and Hansen 2016; Bonhomme and Manresa 2015; Lewis et al. 2023) have been proposed. These methods, which alternate iteratively between estimation and clustering steps, offer a promising approach to refine our understanding of these group-specific relationships and contribute to more accurate and robust analyses in empirical research. The estimation of clustered panel data models is, therefore, important for several reasons. First, it allows for the potential presence of unobserved heterogeneity within the dataset, revealing distinct behaviours that are obscured by traditional panel data analyses. By allowing the estimation of group-specific parameters, these approaches facilitate a nuanced understanding of the underlying dynamics of socio-economic phenomena. Secondly, these methods are particularly valuable when clusters are unknown a priori, as they avoid imposing predetermined grouping structures on the data.

In this paper, we propose a new approach with the aim of extending the clusterwise regression approaches for panel data of Bonhomme and Manresa (2015) and Lewis et al. (2023) to account for spatial dimension. Recognising and incorporating the spatial dimension in clustered panel data models becomes extremely important for a more accurate representation of the relationships between variables. In particular, following previous studies on spatio-temporal clustering (e.g. see Pham 2001;

<figure><img src="image_5.png" /></figure>
Coppi et al. 2010 among the others), we modify the objective function of the afore-mentioned methods by including a spatial penalty. In doing so, the obtained clusters reflect not only the similarity of the relationships between dependent variable and the covariates, but also the similarity structure due to unobserved common characteristics implied by the spatial proximity. By incorporating spatial dimension into clustering, we propose a new procedure that makes clusterwise regression methods for panel data consistent with the fundamental principle of the first law of geophysics, which states that everything is related to everything else, but near things are more related than distant things (Tobler 1970). Moreover, as in Lewis et al. (2023), we include fuzziness in the problem to account for uncertainty in the clustering structure.

We notice that spatial interactions are commonly considered in spatial panel data models. Common types of spatial interactions are the inclusion of a spatially lagged dependent variable, spatial lagged error terms or lagged covariates. Elhorst (2003) proposed an ML estimator of the spatial lag and spatial error models, while Lee Lf (2010) proposed a simple bias correction procedure to improve estimates for this class of models. Yu et al. (2008) constructed a bias corrected estimator for a dynamic model. Another approach for addressing spatial dimension is by spatial fixed effects, that are used to deal with spatially correlated omitted variables. This approach is effective for dealing with group-wise spatial dependence (e.g. see Anselin and Arribas-Bel 2013). Nevertheless, all these approaches assume, as in the classical regression analysis, that all the spatial units have the same relationship between the dependent variable and the covariates. Aquaro et al. (2021) proposed an approach for estimating spatial heterogeneous coefficients but under the assumption that the cluster structure is known. On the contrary, clusterwise regression techniques assume that the cluster structure is unknown and estimate it along with cluster-specific parameters.

In the context of clusterwise regression, the spatial dimension has recently been studied by different papers (e.g. Billé et al. 2017; Lee et al. 2017; Sugasawa and Murakami 2021) but mainly in the case of cross-sectional data. While most of the previous papers proposed two-step procedures, Sugasawa and Murakami (2021) carried out regression and spatial clustering simultaneously. In particular, Sugasawa and Murakami (2021) extended the mixture regression proposed in Sugasawa (2021) to include spatial clustering, by introducing a Potts (1952)-type spatial penalty. We note, however, that only a few clusterwise regression methods for panel data have considered spatial clustering. In particular, Pan et al. (2020) proposed a Markov random field constrained product partition model to estimate cluster-specific panel data regressions under the assumption of unknown group membership. Although this approach enhances spatial clustering for panel data, different from the group fixed-effects methods proposed in Bonhomme and Manresa (2015) and Lewis et al. (2023), it adopts a parametric approach by assuming a Gaussian likelihood function. Moreover, the method proposed in Pan et al. (2020) assumes binary membership and does not account for fuzziness in the cluster structure. Fuzziness is useful for dealing with and quantifying uncertainty in the clustering process. Moreover, it can be seen a generalisation of standard hard clustering methods, where the membership degrees are binary (Bezdek et al. 1984).
We remark that alternative approaches exist for estimating heterogeneous model coefficients. In particular, in the Geographically Weighted Regression (GWR) and the Spatially Varying Coefficient (SVC) methods, the regression coefficients vary continuously over space (e.g. see Yu et al. 2021; Zhu et al. 2014). For these models, unlike the clusterwise regression approaches, there is no need to define the number of clusters, and the parameters vary continuously over space. Nevertheless, we note that this flexibility does not come without costs. In particular, since methods produce a unique set of coefficients for each spatial unit, the results can be harder to summarise and interpret compared to the clear-cut clusters from clusterwise regression-based approaches. In general, the use of clusterwise regression methods is better suited when the focus of the analysis is on the understanding of more aggregate regional patterns rather than highly localised variation, or when researchers do expect that coefficients of neighbouring locations are very similar.

In summary, this paper contributes to two different strands of the literature. First, considering the panel data literature for group fixed-effect estimation (Bonhomme and Manresa 2015; Lewis et al. 2023), we propose a more general approach that leverages both fuzziness and spatial clustering. The proposed new method can be of relevance to researchers interested in using fuzzy group fixed-effect methods (as in Lewis et al. 2023), but want to leverage spatial dimension for clustering units. Indeed, as demonstrated in Anselin and Arribas-Bel (2013), the adoption of cluster-specific fixed effects is particularly useful for dealing with spatial dependencies of unknown form if spatially contiguous units belong to the same cluster. Second, we contribute to the literature of spatio-temporal clustering, in that we extend the spatio-temporal fuzzy clustering algorithm of Coppi et al. (2010) to the case of panel data regression.

To show the effectiveness of the proposed approach, two empirical applications are presented, using different spatial datasets from the literature. The first dataset concerns cigarette consumption in the US states (Baltagi and Levin 1992), while the second concerns non-life insurance consumption in Italian provinces (Millo and Carmeci 2011). These applications illustrate the versatility and applicability of the proposed methodology in revealing spatial patterns and refining our understanding of economic dynamics in different contexts.

The paper is structured as follows. Section 2 discusses in more detail the group fixed-effects estimators proposed by previous studies, which are the starting point of our proposal. Section 3 presents the proposal of this paper, i.e. a group fixed-effects estimator with spatial clustering and the algorithmic implementation. Sections 4.1 and 4.2 discuss two different applications to real datasets, on cigarette consumption in the US and on the demand for insurance products in the Italian provinces, respectively. Section 4.3 shows the results obtained with the benchmark methods, while Sect. 5 concludes with concluding remarks and future research directions.
## 2 Overview of the existing group fixed-effects models

Bonhomme and Manresa (2015) proposed the grouped fixed-effects estimator, which is a form of k-means regression for panel data, to study grouped patterns of heterogeneity. Specifically, given N ($i = 1, ..., N$) statistical units observed for T ($t = 1, ..., T$) times, the authors propose the following linear model

$$y_{it} = \sum_{g=1}^{G} \gamma_{ig} x'_{it} \theta_{g} + v_{it}, \quad (1)$$

where $y_{it}$ is the value of the dependent variable for the i-th unit at time t, and $x_{it}$ is the K-dimensional vector of the independent variables for the i-th unit at time t, which can include both strictly exogenous regressors and lagged outcomes and $v_{it}$ the error term. Moreover, $\gamma$ is the binary membership matrix with generic element $\gamma_{ig}$ related to the membership of the unit i to group g and $\theta_g$ the set of K regression coefficients group g. We assume that the elements of the matrices $\mathbf{x} = (x_{it} : i = 1, ..., N, t = i = 1, ..., T)$ and $\mathbf{v} = (v_{it} : i = 1, ..., N, t = i = 1, ..., T)$ are contemporaneously uncorrelated and independently distributed across the i's. Furthermore, group membership variables $\gamma$'s are unknown and need to be estimated.

The grouped fixed-effects estimator of Bonhomme and Manresa (2015) is defined as follows:

$$f(\theta, \gamma) = \sum_{g=1}^{G} \sum_{i=1}^{N} \gamma_{ig} \|y_i - x'_i \theta_g\|^2, \quad (2)$$

where $y_i = (y_{i1}, ..., y_{iT})$ and $x_i = (x_{i1}, ..., x_{iT})$ are T-dimensional vectors and

$$\|y_i - x'_i \theta_g\|^2 = \sum_{t=1}^{T} (y_{it} - x'_{it} \theta_g)^2.$$

The estimation procedure is attained by seeking for $\hat{\theta}$ and $\hat{\gamma}$ solving a minimization problem based on f in (2), i.e.

$$f(\hat{\theta}, \hat{\gamma}) = \min_{(\theta, \gamma)} f(\theta, \gamma). \quad (3)$$

under the constraints

$$\sum_{g=1}^{G} \gamma_{ig} = 1, \quad \gamma_{ig} \ge 0. \quad (4)$$

For given values of $\theta_g$, the optimal group assignment $g^*$ for each individual unit leads to

$$\hat{\gamma}_{ig}(\hat{\theta}_g) = \gamma_{ig^*}(\theta_{g^*}) = \underset{g \in \{1, ..., G\}}{\operatorname{argmin}} \|y_i - x'_i \theta_g\|^2, \quad (5)$$
which means that we assign the i-th unit to g*-th cluster minimising the squared
Euclidean distance of the residuals, where the g-th group prototype is given by the
predicted value using the parameters of the g-th group $x'_{it}\theta_g$.

Given the grouped structure, the resulting estimator of $\theta$ can then be written as

$$ \hat{\theta}_g(\hat{\gamma}_{ig}) = \theta_{g^*}(\gamma_{ig^*}) = \underset{\theta_g}{\operatorname{argmin}} \sum_{i=1}^{N} \|y_i - x'_i \theta_g\|^2, \quad (6) $$

that is the classical Ordinary Least Square (OLS) estimator within each group. We notice that Bester and Hansen (2016) studied a similar problem but with the assumption of known membership. According to Bester and Hansen (2016), group effects estimators may be thought of as intermediate to pooled and fixed-effects estimators, and allow for improved precision under the assumption that there is no misspecification in the structure of heterogeneity.

The algorithm proposed in Bonhomme and Manresa (2015) alternates between two steps; in the first one (5), each i-th individual unit is assigned to the group g, while in the second step (6), an updated estimate of the parameters $\theta$ is computed using OLS. The two steps are alternate until convergence, that is the objective function (3) of two sequential iterations does not reduce much.

To account for time-invariant fixed effects, the model (1) can be modified as

$$ y_{it} = \sum_{g=1}^{G} \gamma_{ig} x'_{it} \theta_g + \eta_i + v_{it} \quad (7) $$

where $\eta_i$ is the unit-specific fixed effect associated with i-th unit. The model can be simply estimated using demeaned variables, that is

$$ (y_{it} - \bar{y}_i) = \sum_{g=1}^{G} \gamma_{ig} (x'_{it} - \bar{x}_i) \theta_g + (v_{it} - \bar{v}_i). \quad (8) $$

Parameters in this setting can still be estimated with OLS, so the iterative algorithm does not change compared to (1).

However, real-world data exhibit complex patterns, so it is reasonable to assume that the same statistical unit may belong to multiple clusters simultaneously. On the one hand, hard clustering simplifies this complexity by enforcing strict exclusivity, potentially overlooking instances with ambiguous cluster memberships. On the other hand, fuzzy clustering accommodates this ambiguity by assigning membership values that reflect the degree of belonging to each cluster. In this regard, Lewis et al. (2023) proposed a fuzzy extension of the Bonhomme and Manresa (2015) group fixed-effect estimator. The fuzzy group fixed-effect objective function introduces a regularisation to the estimated group membership function so that it is no longer binary. The function f in (2) becomes

$$ f_m(\theta, \gamma) = \sum_{g=1}^{G} \sum_{i=1}^{N} \gamma_{ig}^m \|y_i - x'_i \theta_g\|^2, \quad (9) $$
where $m > 1$ is the fuzziness parameter, while $\gamma_{ig}$ represents the membership degree of the $i$-th unit to the $g$-th group.

The optimization problem in (3) is rewritten as follows

$$f_m(\hat{\theta}, \hat{\gamma}) = \min_{(\theta, \gamma)} f_m(\theta, \gamma), \qquad (10)$$

under the constraints given in (4).

A well-known result in fuzzy clustering is that as $m \to 1^+$, the clustering problem converges to hard clustering (Bezdek et al. 1984). In the context of cross-sectional clustered regression, this result has been shown by Hathaway and Bezdek (1993), while Lewis et al. (2023) demonstrated the same result in the case of fuzzy group fixed-effect regression. Lewis et al. (2023) highlight that the optimal membership degree solving the objective function (10) is given by

$$\hat{\gamma}_{ig} = \sum_{h=1}^{G} \left[ \frac{\|y_i - x'_i \theta_g\|}{\|y_i - x'_i \theta_h\|} \right]^{\frac{2}{m-1}}. \qquad (11)$$

Parameter estimates are obtained through Weighted Least Squares (Hathaway and Bezdek 1993), where the weight of each $i$-th unit for the $g$-th cluster regression is given by the membership degree $\gamma_{ig}$. Indeed, the WLS estimates equal to

$$\hat{\theta}_g = \left( \sum_{i=1}^{N} \gamma_{ig}^m \sum_{t=1}^{T} x'_{it} x_{it} \right)^{-1} \sum_{i=1}^{N} \gamma_{ig}^m \sum_{t=1}^{T} x'_{it} y_{it}, \qquad (12)$$

where we assume that the membership degree is constant over time. Notice that for binary membership, the WLS (12) coincides with group-wise OLS.

Fixed effects like the $\eta$'s in (7) can be easily included also in the fuzzy specification. Therefore, the fuzzy group fixed-effect estimator can be solved iteratively with a similar algorithm to the hard version proposed in Bonhomme and Manresa (2015). The main difference is in the group assignment step, which is replaced by the membership degree updating step in the case of fuzzy group fixed-effect regression.

### 3 Group fixed-effects estimation with spatial clustering

Previous approaches overlook the spatial dimension, neglecting the inherent relationships and dependencies that exist between adjacent units. This oversight can lead to suboptimal results, particularly in applications where spatial factors play a relevant role. In such cases, the presence of fixed effects common to groups of spatially neighbouring units is of particular importance, as they capture common characteristics between spatially neighbouring units and thus reflect the interdependence that arises from geographical proximity or spatial relationships (e.g. Anselin and Arribas-Bel 2013; Millo and Carmeci 2011; Von Lyncker and Thoennessen 2017). Neglecting these fixed effects can lead to biased estimates and incomplete representations of the underlying structure in the data.
One way to extend the model in (10) to spatial clustering is to add an appropriate spatial penalty term in the minimisation problem. In the following, we say that two units are adjacent if and only if they share a common boundary. The spatial penalty captures spatial autocorrelation and ensures that contiguous spatial units are grouped appropriately. For this purpose, we introduce a square symmetric matrix of order N, denoted P. The diagonal elements of P are set to zero, while the generic off-diagonal element $p_{ij}$, for units i, j where i ≠ j, is given by

$$p_{ij} = \begin{cases} 1 & \text{if unit } i \text{ is contiguous to unit } j, \\ 0 & \text{otherwise.} \end{cases}$$

We notice that, although the binary contiguity matrix is the most commonly adopted, several possible matrices describing the spatial information can be built (Fouedjio 2016; Mattera 2022) without loss of generality. A suitable spatial penalty term proposed by the previous studies (Pham 2001; Coppi et al. 2010; D'Urso et al. 2023) is

$$\mathrm{Pen}_m(\gamma, \mathbf{P}) = \frac{\beta}{2} \sum_{i=1}^{N} \sum_{g=1}^{G} \gamma_{ig}^m \sum_{j=1}^{N} p_{ij} \sum_{h \in G_g} \gamma_{jh}^m, \quad (13)$$

where $\beta \ge 0$ is the penalty coefficient, which regulates the impact of the spatial penalization term, and $G_g = \{1, ..., g-1, g+1, ..., G\}$. The spatial penalty term in (13) is based on the following assumption. If a spatial unit i belongs to cluster g with a high membership degree, the (13) forces the neighbouring spatial units to have high membership degrees in cluster g. Therefore, the penalty (13) determines a spatially smoothed membership degree, assuming that approximately similar features often characterise neighbouring areas. Indeed, spatial proximity has been used in the previous studies as a way of instrumenting for similar unobserved characteristics (Mattera and Franses 2023). It is important to note that the inclusion of a spatial penalty term does not necessarily force the resulting clusters to consist only of adjacent spatial units, but it does tend to encourage the formation of clusters with contiguous areas. Indeed, the spatial penalty term does not dictate that adjacent areas must necessarily belong to the same cluster; rather, it seeks to identify clusters that include adjacent areas. However, in situations where a few spatial units have characteristics that are markedly different from those of their neighbours, the clustering process does not prevent these specific units from being assigned to different clusters.

Following Coppi et al. (2010) and D'Urso et al. (2019), we incorporate the spatial penalty in the problem (10) by means of the following optimization problem

$$\min_{(\theta, \gamma)} [f_m(\theta, \gamma) + \mathrm{Pen}_m(\gamma, \mathbf{P})], \quad (14)$$

under the constraints in (4) and where $f_m$ and $\mathrm{Pen}_m$ are defined as in (9) and (13), respectively.

**Remark** We remark that the proposed spatial fuzzy GFE encompasses the methods proposed in Bonhomme and Manresa (2015) and Lewis et al. (2023). As shown
in Fig. 1, by setting $\beta = 0$, we remove the spatial penalty and obtain the Lewis et al. (2023) model, while further constraining $m = 1$, we get the Bonhomme and Manresa (2015) approach.

It is worth noting that for the choice of $m$ and $G$, we can consider appropriate cluster validity criteria. Similarly, for the optimal choice of $\beta$, we can explore appropriate methods. For each spatial unit $i$ and each cluster $g$, the objective is to minimise the sum of the membership degrees of contiguous spatial units (as specified in the matrix $P$) in all clusters except cluster $g$ (represented by $G_g = \{1, \dots, g-1, g+1, \dots, G\}$). This constraint ensures the spatial homogeneity of the clusters while maximising the internal cohesion based on the feature vectors $x_i$. The parameter $\beta$ plays a crucial role in balancing the trade-off between encouraging contiguous/neighbouring spatial units to belong to the same cluster (for $\beta \gg 0$) and neglecting spatial regularisation altogether (for $\beta = 0$). The optimal solution to the minimization problem (14) is given by Proposition 1.

**Proposition 1** The optimal solution of the problem (14) with constraints in (4) is the couple $(\hat{\theta}, \hat{\gamma})$, such that

$$ \hat{\theta}_g = \left( \sum_{i=1}^{N} \gamma_{ig}^m \sum_{t=1}^{T} x'_{it} x_{it} \right)^{-1} \sum_{i=1}^{N} \gamma_{ig}^m \sum_{t=1}^{T} x'_{it} y_{it}, \quad (15) $$

is the regression coefficients of the $g$-th component of $\hat{\theta}$ and

$$ \hat{\gamma}_{ig} = \left[ \sum_{h=1}^{G} \frac{\|y_i - x'_i \theta_g + \beta \sum_{j=1}^{N} \sum_{h \in G_g} p_{ij} \gamma_{jh}^m \|}{\|y_i - x'_i \theta_h + \beta \sum_{j=1}^{N} \sum_{l \in G_l} p_{ij} \gamma_{jl}^m \|} \right]^{\frac{2}{m-1}} \quad (16) $$

is the component of $\hat{\gamma}$ describing the membership degree of unit $i$ for group $g$.

**Proof** The Lagrangian of the considered constrained optimization problem is

$$ \mathcal{L} = \sum_{g=1}^{G} \sum_{i=1}^{N} \gamma_{ig}^{m} \left[ \|y_i - x'_i \theta_g\|^{2} + \frac{\beta}{2} \sum_{j=1}^{N} \sum_{h \in G_g} p_{ij} \gamma_{jh}^{m} \right] + \sum_{i=1}^{N} \lambda_i \left( 1 - \sum_{g=1}^{G} \gamma_{ig} \right). \quad (17) $$

<figure><img src="image_6.png" /><figcaption>**Fig. 1** Special cases of the proposed approach. In writing $A \to B$, we indicate that $B$ is a subcase of $A$. In particular, by setting $\beta = 0$ we remove the spatial penalty from the proposed general model and obtain the Lewis et al. (2023) model, while further constraining $m = 1$, we get the Bonhomme and Manresa (2015) approach</figcaption></figure>
To find the optimal regression parameters, we need to take the first derivative of the Lagrangian with respect to $\theta_g$ and set it equal to zero. However, we recognise that the solution to this problem is identical to the one obtained for (10), for which we already know that $\hat{\theta}$ is as in (12), i.e. as the well-known WLS estimator in (15).

To find the optimal membership degree, we take the partial derivative of (17) to $\gamma_{ig}$, that is

$$ \frac{\partial \mathcal{L}}{\partial \gamma_{ig}} = m \gamma_{ig}^{m-1} \left[ \|y_i - x'_i \theta_g\|^2 + \beta \sum_{j=1}^{N} \sum_{h \in G_g} p_{ij} \gamma_{jh}^m \right] - \lambda_i = 0, $$

which yields to

$$ \gamma_{ig} = \left[ \frac{m \|y_i - x'_i \theta_g\|^2 + \beta \sum_{j=1}^{N} \sum_{h \in G_g} p_{ij} \gamma_{jh}^m}{\lambda_i} \right]^{-1/(m-1)}. \quad (18) $$

By employing the constraint in (14), we obtain

$$ \sum_{g=1}^{G} \gamma_{ig} = \sum_{g=1}^{G} \left[ \frac{m \|y_i - x'_i \theta_g\|^2 + \beta \sum_{j=1}^{N} \sum_{h \in G_g} p_{ij} \gamma_{jh}^m}{\lambda_i} \right]^{-1/(m-1)} = 1 $$

thus, solving for $\lambda_j$, we get

$$ \lambda_i^{-1/(m-1)} = \sum_{g=1}^{G} \left[ m \|y_i - x'_i \theta_g\|^2 + \beta \sum_{j=1}^{N} \sum_{h \in G_g} p_{ij} \gamma_{jh}^m \right]^{-1/(m-1)}, \quad (19) $$

since $\lambda_i$ does not depend on $g$, it can be taken out from the summation. By substituting (19) in (18), we obtain (16). □

We note, however, that the obtained closed-form solution to the problem has a complex shape, since we need the membership degree $\gamma_{ig}$ to estimate $\theta_g$ with WLS, but we need $\theta_g$ to estimate $\gamma_{ig}$. The solution has to be iterative. Following Coppi et al. (2010) and D’Urso et al. (2019), and in line Bonhomme and Manresa (2015), we propose an algorithm for estimating the fuzzy grouped fixed-effect model with spatial clustering (spatial fuzzy GFE), which works as follows: given the choice of the tuning parameters $m$, $G$ and $\beta$, we consider an initial clustering (e.g. based on spatial coordinates, see Sugasawa and Murakami 2021). Then, the algorithm alternates the estimation and clustering assignment steps until convergence. Note that the iterative algorithm converges at least to a local optimum. To reduce the risk of hitting local optima, more than one random start is recommended. Regarding the theoretical properties of the proposed fuzzy GFE estimator, we note that the introduction of spatial clustering does not modify the properties of the classical GFE estimator studied in Bonhomme and Manresa (2015) and Lewis et al. (2023).
To choose $\beta$ and $G$, following Coppi et al. (2010) and D'Urso et al. (2022), we suggest to maximise a measure of the spatial autocorrelation of the clusters, taking into account the fuzzy membership matrix $\gamma$ and the spatial proximity matrix $\mathbf{P}$. Let $\hat{\mathbf{E}}$ be the $(N \times G)$ dimensional matrix containing the squared Euclidean distance of each $i$-th unit to the $g$-th cluster prototype $\hat{e}_{ig} = \|y_i - x'_i\hat{\theta}_g\|^2$, the spatial autocorrelation measure for the $g$-th cluster is given by

$$ \rho_g = \frac{\mathrm{tr}\left[\hat{\mathbf{E}}'_g \gamma_g^{\frac{1}{2}} \mathbf{P}_g^{\frac{1}{2}} \gamma_g^{\frac{1}{2}} \hat{\mathbf{E}}_g\right]}{\mathrm{tr}\left[\hat{\mathbf{E}}_g \gamma_g^{\frac{1}{2}} \mathrm{diag}(\mathbf{PP}) \gamma_g^{\frac{1}{2}} \hat{\mathbf{E}}_g\right]}, \quad (20) $$

where $\gamma_g$ is the order $N$ square diagonal matrix of the membership degrees of cluster $g$, and $\mathbf{P}$ is the spatial contiguity matrix. Similar to the fuzzy Moran index (D'Urso et al. 2019), $\rho_g$ ranges between -1 and 1. A value of 1 indicates perfect positive spatial autocorrelation, i.e. neighbouring units have similar values, 0 indicates no autocorrelation, i.e. units are spatially randomly distributed, and -1 indicates perfect negative spatial autocorrelation, i.e. neighbouring units have dissimilar values. Thus, the higher the value of these statistics, the better the geographical allocation of units to clusters. To determine the overall measure of spatial autocorrelation for the partition obtained, we can calculate the weighted mean of the measure, i.e.

$$ \rho = \frac{\sum_{g=1}^{G} \rho_g \sum_{i=1}^{N} \gamma_{ig}}{N}. \quad (21) $$

Following Coppi et al. (2010), after running the clustering models for increasing values of spatial penalty $\beta$, we choose the value maximising overall spatial autocorrelation $\rho$ in (21).

# 4 Empirical analysis to real data

## 4.1 Cigarette consumption in the US states

Studying the dynamics of cigarette consumption in response to taxation policies provides valuable insights into the effectiveness of price-based interventions in determining smoking initiation and promoting cessation efforts. In this regard, Baltagi and Levin (1986) has shown some of the pitfalls of previous studies that relied on either time-series regression of a given state or cross-sectional regression for a given year, and introduced the “bootlegging” effect into the modelling equation. Bootlegging is explicitly modelled using the “neighbourhood” price and has also been found to be significant in explaining cigarette consumption in the US. Baltagi and Levin (1992) updated the Baltagi and Levin (1986) results and confirmed the existence of positive negative price effects on consumption, suggesting that taxation can be an effective tool for reducing smoking behaviour in the population.
Therefore, the effectiveness of cigarette taxation as a policy instrument depends on a reliable estimate of the price elasticity of demand for cigarettes. However, we argue that such elasticity may differ across geographic states due to various socio-economic, cultural and regulatory factors. For example, states with higher income levels and stronger anti-smoking campaigns may have greater price sensitivity among consumers, leading to a higher price elasticity of demand. Conversely, countries with lower-income levels or cultural norms that are more tolerant of smoking may have lower price elasticity, indicating less responsiveness to changes in cigarette prices. In addition, variations in tobacco control policies, such as differences in tobacco advertising regulations or smoking cessation programmes, may also influence the price elasticity of demand across subgroups of the US states. This source of heterogeneity cannot be captured by standard approaches, even those that control for spatial effects. Conversely, the method proposed in Sect. 3 can highlight the existence of such patterns, which may be useful for policy perspectives.

The sample consists of a panel dataset from 46 American states over the period 1963 to 1988, which is the same dataset studied in Baltagi and Levin (1992). The study areas are shown in Fig. 2. Oregon and Colorado were excluded due to missing data.

Baltagi and Levin (1986) and Baltagi and Levin (1992) estimate the following econometric model for cigarette demand in the US

$$ \ln(C_{it}) = \theta_0 + \theta_1 \ln(C_{it-1}) + \theta_2 \ln(P_{it}) + \theta_3 \ln(Y_{it}) + \theta_4 \ln(Pn_{it}) + \varepsilon_{it}, \quad (22) $$

where $C_{it}$ is real per capita sales of cigarettes by persons of smoking age (14 years and older), measured in packs of cigarettes per head, $P_{it}$ is the average retail price of

US states under study

<figure><img src="image_7.png" /><figcaption>Fig. 2 Selected study area</figcaption></figure>
a pack of cigarettes measured in real terms, $Y_{it}$ is real per capita disposable income and $Pn_{it}$ denotes the minimum real price of cigarettes in any neighbouring state.

Before applying the spatial fuzzy GFE approach, we need to choose the fuzzi-
fier parameter *m*, the number of clusters *G* and the value of the spatial penalty *β*.
About the fuzzifier parameter *m*, we fix *m* = 2. There are several reasons behind
this choice. The previous studies suggest that setting *m* = 2 in fuzzy *k*-means
clustering strikes a balance between interpretability, computational efficiency
and effective cluster separation, making this a widely adopted choice in practice
(D'Urso et al. 2015; Pal and Bezdek 1995; Wang et al. 2007). Indeed, setting
*m* = 2 ensures that the clusters maintain a balance between soft assignments and
clear separations. Lowering *m* artificially reduces fuzziness and tends to move
the model closer to hard clustering, where each data point belongs more defini-
tively to one cluster. While this may lead to sharper separation between clusters,
it also limits the model's ability to capture uncertainty and gradual transitions
between regions where cluster boundaries are not well-defined. Moreover, the
choice of *m* = 2 is consistent with the squared error loss used for least squares
estimation. The quadratic specification ensures that larger distances are penal-
ised more heavily, meaning that data points far from the cluster prototype have
a reduced influence on the estimation process. By assigning stronger penalties to
larger distances, the model ensures that the clustering is driven by points that are
more representative of the local structure, while still accounting for the presence
of outliers in a balanced way.

After fixing *m*, we optimally choose both *G* and *β*. For this aim, as discussed in Sect. 3, we follow a strategy similar to Coppi et al. (2010) based on maximising the intra-cluster spatial autocorrelation (21). Considering a grid of *G* = 2, 3, ..., 10 and *β* = 0, 0.1, ..., 2, we choose the best combination of tuning parameters that maximise the overall spatial autocorrelation. The best choice is *G* = 2 and *β* = 0.9. Figure 3 shows the values of *β* for different values of *G*. We report the cases *G* = 1, ..., 6 for a better visualisation. The curves with different colours refer to the *ρ* value obtained with different numbers of clusters *G*. In summary, we choose *β* = 0.9 because the spatial autocorrelation of the clusters does not increase after this cut-off point, and this result holds for all the considered values of *G*. Interestingly, the largest value of *ρ* is obtained with *G* = 2 but the difference with other numbers of clusters is negligible. Thus, we choose *G* = 2 by maximising *ρ* but ensuring at the same time parsimony.

The clustering structure found by the algorithm is shown in Fig. 4, while Table 1
shows the regression estimates, and Table 2 shows the membership degree matrix.
The first column of Table 1, labelled “No clustering”, shows the results of the within
OLS regression with fixed effects but without clustering effects, while the other two
columns give the estimate for the two clusters obtained with the spatial fuzzy GFE
model.

Figure 4 also shows fuzzy states, i.e. those with a degree of cluster membership lower than 0.7 (see also Table 2). Indeed, the chosen cut-off point of 0.7 for a partition into two clusters for the membership degrees is compatible with the indications suggested in the literature (Dembele and Kastner 2003). In summary, the spatial fuzzy GFE highlights the presence of two main clusters, namely the eastern US
Selecting spatial penalty

<figure><img src="image_8.png" /><figcaption>**Fig. 3** Values of *ρ* for different values of the spatial penalty *β* and clusters *G*. Lines with different colours refer to alternative values of *G*. Cigarette consumption dataset</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 4** Clustering structure found by the spatial fuzzy GFE method. Blue states are denoted as fuzzy, that is the US states with membership degree to the clusters lower than 0.7</figcaption></figure>

<figure><img src="image_10.png" /></figure>
states and the western states. Some states in the centre of the US are considered fuzzy, i.e. not belonging to either group.

Overall, the regression estimates obtained by the spatial fuzzy GFE are consist-
ent with the results of the previous studies. Indeed, we find that lagged consump-
tion values θ₁ have a positive effect in explaining current cigarette consumption, thus
documenting the smoking habit component. Moreover, the price elasticity given by
the coefficient θ₂ is negative and statistically significant, and also the magnitude of
the estimates is very close between the clusters as well as compared to the estimates
without clustering. Instead, a first difference can be found in the income elasticity,
i.e. the coefficient θ₃, which is positive in both groups but much larger for cluster 2
than for cluster 1. Indeed, it can be noted that the θ₃ estimated without clustering
effects is close to an average of the coefficients estimated in both clusters, while this
difference in magnitude may be relevant from a policy perspective. Finally, the most
significant difference can be found in the coefficient θ₄, that is the smuggling effect.
This effect is significant according to the within OLS estimates, but the spatial fuzzy
GFE highlights that this effect is not statistically significant for units placed in clus-
ter 2. Thus, we find new evidence that the bootlegging effect is positive and signifi-
cant only for some subgroups of units and not for all US states.

One possible explanation for this finding is that the US states tend to be
smaller and more densely populated than their western counterparts. This means
that residents of eastern states may have easier access to neighbouring states
where cigarette prices may be lower due to differences in tax policies. In con-
trast, western US states often have large distances between them, making it less
practical for residents to travel across state lines to buy cheaper cigarettes. In this

<figure><img src="image_11.png" /><figcaption>**Fig. 5** Median real income in the US states for the considered time period. The figure shows that western states are less densely populated compared to eastern states, especially those in the north-east</figcaption></figure>
US states: Real Income

<figure><img src="image_12.png" /><figcaption>**Fig. 6** Median real income in the US states for the considered time period. The figure shows that, with some exceptions, eastern states are wealthier compared to those in the west</figcaption></figure>

<table><caption>**Table 1** Spatial fuzzy GFE with $G = 2$, $\beta = 0.9$ and $m = 2$</caption><thead><tr><th></th><th>No clustering</th><th>Cluster 1</th><th>Cluster 2</th></tr></thead><tbody><tr><td>θ<sub>0</sub></td><td>0.17<sup>***</sup><br>(0.03)</td><td>0.26<sup>***</sup><br>(0.03)</td><td>0.14<sup>***</sup><br>(0.03)</td></tr><tr><td>θ<sub>1</sub></td><td>0.96<sup>***</sup><br>(0.01)</td><td>0.95<sup>***</sup><br>(0.01)</td><td>0.94<sup>***</sup><br>(0.01)</td></tr><tr><td>θ<sub>2</sub></td><td>-0.10<sup>***</sup><br>(0.01)</td><td>-0.11<sup>***</sup><br>(0.01)</td><td>-0.10<sup>***</sup><br>(0.02)</td></tr><tr><td>θ<sub>3</sub></td><td>0.03<sup>***</sup><br>(0.01)</td><td>0.02<sup>***</sup><br>(0.01)</td><td>0.06<sup>***</sup><br>(0.01)</td></tr><tr><td>θ<sub>4</sub></td><td>0.04<sup>***</sup><br>(0.01)</td><td>0.06<sup>***</sup><br>(0.01)</td><td>0.00<br>(0.02)</td></tr><tr><td>Adj. R<sup>2</sup></td><td>0.97</td><td>0.97</td><td>0.97</td></tr></tbody></table>

Standard errors are in parenthesis

***p < 0.001; **p < 0.01; *p < 0.05

regard, Fig. 5 shows the spatial heterogeneity of the median population density in the US states. The figure highlights that western states are less densely populated compared to eastern states, especially those in the north-east. The states in the Midwest, on the other hand, are those with lower levels of population density. In the west, the main exception is the state of California. In addition, eastern US states have historically been more interconnected (Nelson and Rae 2016), with
<table><caption>**Table 2** Membership degree matrix resulting from the spatial fuzzy GFE</caption><thead><tr><th>US states</th><th>Cluster 1</th><th>Cluster 2</th></tr></thead><tbody><tr><td>Alabama</td><td>0.0160</td><td>0.9840</td></tr><tr><td>Arizona</td><td>0.9879</td><td>0.0121</td></tr><tr><td>Arkansas</td><td>0.4386</td><td>0.5614</td></tr><tr><td>California</td><td>0.9813</td><td>0.0187</td></tr><tr><td>Connecticut</td><td>0.0123</td><td>0.9877</td></tr><tr><td>Delaware</td><td>0.0099</td><td>0.9901</td></tr><tr><td>Florida</td><td>0.0128</td><td>0.9872</td></tr><tr><td>Georgia</td><td>0.0056</td><td>0.9944</td></tr><tr><td>Idaho</td><td>0.9885</td><td>0.0115</td></tr><tr><td>Illinois</td><td>0.3573</td><td>0.6427</td></tr><tr><td>Indiana</td><td>0.0522</td><td>0.9478</td></tr><tr><td>Iowa</td><td>0.7929</td><td>0.2071</td></tr><tr><td>Kansas</td><td>0.8766</td><td>0.1234</td></tr><tr><td>Kentucky</td><td>0.0826</td><td>0.9174</td></tr><tr><td>Louisiana</td><td>0.4040</td><td>0.5960</td></tr><tr><td>Maine</td><td>0.0242</td><td>0.9758</td></tr><tr><td>Maryland</td><td>0.0051</td><td>0.9949</td></tr><tr><td>Massachusetts</td><td>0.0035</td><td>0.9965</td></tr><tr><td>Michigan</td><td>0.1248</td><td>0.8752</td></tr><tr><td>Minnesota</td><td>0.9019</td><td>0.0981</td></tr><tr><td>Mississippi</td><td>0.1368</td><td>0.8632</td></tr><tr><td>Missouri</td><td>0.5631</td><td>0.4369</td></tr><tr><td>Montana</td><td>0.9921</td><td>0.0079</td></tr><tr><td>Nebraska</td><td>0.9290</td><td>0.0710</td></tr><tr><td>Nevada</td><td>0.9828</td><td>0.0172</td></tr><tr><td>New Hampshire</td><td>0.0373</td><td>0.9627</td></tr><tr><td>New Jersey</td><td>0.0079</td><td>0.9921</td></tr><tr><td>New Mexico</td><td>0.9575</td><td>0.0425</td></tr><tr><td>New York</td><td>0.0037</td><td>0.9963</td></tr><tr><td>North Dakota</td><td>0.9799</td><td>0.0201</td></tr><tr><td>Ohio</td><td>0.0075</td><td>0.9925</td></tr><tr><td>Oklahoma</td><td>0.8168</td><td>0.1832</td></tr><tr><td>Pennsylvania</td><td>0.0058</td><td>0.9942</td></tr><tr><td>Rhode Island</td><td>0.0407</td><td>0.9593</td></tr><tr><td>South Carolina</td><td>0.0256</td><td>0.9744</td></tr><tr><td>South Dakota</td><td>0.9838</td><td>0.0162</td></tr><tr><td>Tennessee</td><td>0.0987</td><td>0.9013</td></tr><tr><td>Texas</td><td>0.7314</td><td>0.2686</td></tr><tr><td>Utah</td><td>0.9865</td><td>0.0135</td></tr><tr><td>Vermont</td><td>0.0182</td><td>0.9818</td></tr><tr><td>Virginia</td><td>0.0085</td><td>0.9915</td></tr><tr><td>Washington</td><td>0.9505</td><td>0.0495</td></tr><tr><td>West Virginia</td><td>0.0061</td><td>0.9939</td></tr><tr><td>Wisconsin</td><td>0.5505</td><td>0.4495</td></tr><tr><td>Wyoming</td><td>0.9872</td><td>0.0128</td></tr></tbody></table>

Cigarette consumption in the US dataset
higher rates of interstate travel. As a result, residents of eastern states may be more inclined to cross state lines to take advantage of lower cigarette prices in neighbouring states. On the other hand, western US states may have less interstate travel, leading to a weaker influence of neighbouring states' prices on cigarette consumption. This explanation is in accordance with the regression results (see Table 1), for which the neighbouring states' prices effect is statistically significant compared to states in the west.

However, also economic disparities between neighbouring states may also play a role. Eastern US states often have greater economic disparities within close proximity. Figure 6 shows the spatial distribution of the median real income between the US states for the considered time period. Interestingly, we find that higher values of income can be found in the states in the east, specifically those in the north-east. The states in the Midwest, on the other hand, are characterised by lower-income levels. Residents in lower-income areas of eastern states may be more sensitive to price differences and thus more likely to seek cheaper cigarettes across state lines. In contrast, economic disparities between neighbouring states in the west may be less pronounced, leading to a reduced impact of prices in neighbouring states on cigarette consumption. Moreover, we highlight that economic differences also correlate with smoking habits. For instance, Lawrence et al. (2007) report that unemployed and workers reporting a low annual household income are more likely to report current, daily and heavy smoking. In this regard, Fig. 6 highlights an interesting difference between west and east in terms of real income. Cultural attitudes towards smoking and tobacco use may differ between eastern and western US states. Eastern states, particularly those with larger urban populations, may have stronger anti-smoking campaigns and regulations, leading to greater price sensitivity among consumers. In contrast, western states, especially those with a history of tobacco production or less stringent tobacco control measures, may have different consumption patterns that are less influenced by prices in neighbouring states. Using data about a similar period of that analysed in this application, Alciati et al. (1998) constructed an indicator of the policy stringency on tobacco usage in the US states. Considering the scores computed for the early 90s, the authors find that states in the Midwest were those with the less stringent policy and that in general states in the east adopted more strict policies against tobacco consumption compared to those in the west. Similar results have been discussed in Botello-Harbaum et al. (2009). Therefore, the distinction between two clusters in Fig. 4 is also consistent with the evidence of the previous studies.

Overall, the combination of geographic proximity, interstate travel patterns, economic disparities and cultural/regulatory factors contribute to the greater importance of neighbouring states' prices on cigarette consumption in the eastern US states compared to those in the west. Such differences are difficult to identify using standard approaches.

<figure><img src="image_13.png" /></figure>
## 4.2 Insurance demand in the Italian provinces

Understanding the factors that drive non-life insurance consumption is crucial for both policymakers and the industry. Policymakers can use this knowledge to design policies that promote greater insurance uptake, while insurance companies can gain insights into consumer behaviour, enabling them to develop tailored marketing strategies and product offerings, improve market penetration and meet customer needs. By identifying these determinants, policymakers and stakeholders could also work together to improve society's overall financial security and resilience.

The previous literature on the determinants of non-life insurance (e.g. see Outreville 1990; Gruber and Poterba 1994) has evolved from the use of cross-sectional analysis to the use of panel data techniques to mitigate the influence of unobservable factors. However, as highlighted by Millo and Carmeci (2011), these regression models remain susceptible to unobserved heterogeneity, mainly due to institutional and geographical factors, which may not be fully captured by simple fixed-effects models. Furthermore, Millo (2016) emphasises that the coefficients of interest may exhibit heterogeneity across spatial units, highlighting the need for nuanced analysis.

The Italian context provides a compelling case study to explore these dynamics. Indeed, on average, the Italian non-life insurance market lags significantly behind its European counterparts. Moreover, an examination of the geography of Italian insurance consumption reveals stark differences between north and south. This discrepancy reflects the broader economic heterogeneity within Italy (Accetturo et al. 2024; Menon et al. 2023; Salvati 2013). Therefore, the study of the heterogeneity of the determinants of non-life insurance consumption in Italy, with an emphasis on the adoption of a spatial perspective, becomes particularly interesting against the background of these pronounced regional disparities.

Building on the findings of Millo and Carmeci (2011), our study examines non-life insurance consumption patterns in 103 Italian provinces from 1998 to 2002, to identify key determinants. Through our methodological approach, we aim to uncover additional nuances and characteristics that may provide further insights through clustering techniques. Compared to the previous application in Sect. 4.1, we now have a short panel with a large cross-sectional dimension. We consider the following model (Millo and Carmeci 2011)

$$ \ln(\text{ppc}_{it}) = \theta_0 + \theta_1 \ln(\text{rgdp}_{it}) + \theta_2 \ln(\text{rbdp}_{it}) + \theta_3 \ln(\text{dens}_{it}) \\ + \theta_4 \text{rirs}_{it} + \theta_5 \ln(\text{agen}_{it}) + \theta_6 \text{school}_{it} + \theta_7 \text{vagr}_{it} \\ + \theta_8 \text{family}_{it} + \theta_9 \text{inef}_{it} + \theta_{10} \text{trust}_{it} + \varepsilon_{it}, \quad (23) $$

where ppc$_{it}$ is the premium per capita. The controls are rgdp$_{it}$ the regional GDP, rbdp$_{it}$ the real bank deposit, dens$_{it}$ is the density of inhabitants, rirs$_{it}$ is the real interest rate on borrowing, agen$_{it}$ is the density of insurance agencies, school$_{it}$ is the share of well-educated people, vagr$_{it}$ share of value added (agricultural), family$_{it}$ is the average number of family members, ineff$_{it}$ is the judicial inefficiency and trust$_{it}$ reflects the share of people declaring to trust others. The selected study area is shown in Fig. 7
<figure><img src="image_14.png" /><figcaption>**Fig. 7** Selected study area: Italian provinces</figcaption></figure>

To apply the spatial fuzzy GFE approach, again fixing the fuzzifier parameter at $m = 2$, we choose the number of clusters $G$ and the value of the spatial penalty $\beta$ by maximising the intra-cluster spatial autocorrelation (21). Considering a grid of $G = 2, 3, ..., 10$ and $\beta = 0, 0.1, ..., 2$, we choose $G = 2$ and $\beta = 0.1$. Figure 8 shows the values of $\beta$ for different choices of $G$. We notice that, for all the values of $G$, the value of the spatial penalty is much lower compared to the empirical application to cigarette consumption shown in Sect. 4.1, and that the value of the total spatial autocorrelation does not increase much from $\beta = 0$ to $\beta = 0.1$, i.e. from $\rho = 0.268$ to $\rho = 0.280$. This result can be easily explained by the fact that even with $\beta = 0$ the spatial GFE identifies clusters of units that are spatially clustered. In other words, the spatial penalty is needed to enhance a natural spatial clustering present in the data, given by the similarity of the parameter estimates in the different clusters. This result is important because it clearly shows that the approach used to select $\beta$ works well, i.e. $\beta$ is close to zero when the spatial penalty is not really needed because a natural spatial clustering already exists. Given that the difference in terms of $\rho$ between the cases $G = 2$ and $G = 5$, we choose the solution with $G = 2$ clusters with the aim of parsimony. In other words, there is no hige difference in terms of $\rho$ between the two partitions with $G = 2$ and $G = 5$, so we find more appropriate the solution in two main clusters. The resulting partition, shown in Fig. 9, distinguishes between the north and the south of Italy, and this is coherent with the economic gap between the two main spatial aggregates of the country (e.g. see Accetturo et al. 2024).
Selecting spatial penalty

<figure><img src="image_15.png" /><figcaption>**Fig. 8** Values of *ρ* for different values of the spatial penalty *β* and clusters *G*. Lines with different colours refer to alternative values of *G*. Insurance consumption dataset</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 9** Clustering structure
find by the spatial fuzzy GFE
method. Cluster 1 is highlighted
in green, while Cluster 2 is in
red</figcaption></figure>
Table 3 shows the regression estimates. The first column of Table 3, labelled “No clustering”, shows the results of the within OLS regression with fixed effects but without clustering effects, while the other two columns give the estimate for the two clusters obtained with the spatial fuzzy GFE model.

Figure 9 perfectly reconstructs the differences between the north and the south of Italy in terms of the determinants of non-life insurance consumption. Indeed, the economic disparities between north and south (Accetturo et al. 2024; Salvati 2013) can affect insurance consumption, as wealthier individuals or regions tend to demand more insurance products. Moreover, insurance demand is partly influenced by cultural attitudes towards risk. From this perspective, it is known that southern Italy is traditionally considered more risk-averse (Guiso and Paiella 2008) due to economic uncertainties, and this could lead to different insurance consumption patterns compared to the north. Furthermore, it is known that more developed financial markets are associated with larger insurance consumption (Li et al. 2007). Therefore, it is reasonable to assume that provinces with more developed financial markets—which are typically found in the north of Italy (A’hearn 2005)—tend to have higher insurance penetration.

The partition in Fig. 9 reflects such differences. The only exceptions in the north are Udine and Gorizia, which have similar determinants to the south, and Sassari, which has similar values to the north. In central Italy, the provinces of Molise and Abruzzo have similar determinants to the southern provinces. Interestingly, no fuzzy units are highlighted in Fig. 9.

Figure 10 shows the differences in terms of distribution between North and South clusters for the dependent variable of the model (23) and for the following covariates: real GDP per capita, population density and schooling. The relevant gaps for all the variables well explain the suitability of the north–south partition shown in Fig. 9.

First, we note that the south regions have a significantly lower demand of non-life insurance products and, as briefly mentioned above, the determinants of this may be found in the profound socio-economic differences in the two parts of the country. Richer regions may have consume more compared to others, and indeed, there is an evident gap between north and south, as documented also in other studies (Accetturo et al. 2024; Salvati 2013). Also, the availability of insurance agencies may play a role in explaining the lower level of insurance demand, and in Fig. 10, we find that the density of insurance is much lower in the South cluster. Finally, a more educated population may better understand the importance of insurance products and demand more. Also in terms of schooling (see Fig. 10c), we find significant differences in the two clusters.

Table 3 highlights the differences in the determinants of non-line insurance. First of all, we notice that some determinants that appear to be significant when considering the full sample analysis (first column of Table 3) disappear when considering the spatial fuzzy GFE estimates. In fact, population density (the coefficient $\theta_3$) seems to be a positive and statistically significant factor for insurance consumption, while clustered results show that such a determinant is significant only for the cluster associated with the northern provinces. Similarly, the trust variable (coefficient $\theta_{10}$) is strongly significant only for the southern

<figure><img src="image_17.png" /></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. 10** Differences in terms of distribution between North and South clusters for the dependent variable of the model (23) and for the following covariates: real GDP per capita, population density and schooling</figcaption></figure>

provinces, while it is not significant for the northern provinces. This result has a particularly strong implication, as the lack of trust in others in the southern provinces explains well the lower consumption in the south compared to the north. This result is not highlighted by standard approaches. Other variables, such as the real interest rate on loans and the density of insurance agencies, are instead statistically significant with the same sign in both clusterings, but with significantly different magnitudes.

## 4.3 Results with benchmark methods

To validate the usefulness of the proposed approach, in what follows we discuss the results obtained with the benchmark methods, that is Bonhomme and Manresa (2015) and Lewis et al. (2023). Both approaches do not employ spatial constraints, while the Lewis et al. (2023) group fixed-effect approach adopts fuzziness, the Bonhomme and Manresa (2015) does not.
<table><caption>**Table 3** Spatial fuzzy GFE with *G* = 2, *β* = 0.1 and *m* = 2. Standard errors are in parenthesis</caption><thead><tr><th></th><th>No clustering</th><th>Cluster 1</th><th>Cluster 2</th></tr></thead><tbody><tr><td>θ<sub>0</sub></td><td>-7.54<sup>***</sup><br>(0.79)</td><td>-5.79<sup>***</sup><br>(0.69)</td><td>-9.91<sup>***</sup><br>(0.86)</td></tr><tr><td>θ<sub>1</sub></td><td>0.91<sup>***</sup><br>(0.08)</td><td>0.91<sup>***</sup><br>(0.07)</td><td>0.98<sup>***</sup><br>(0.09)</td></tr><tr><td>θ<sub>2</sub></td><td>0.43<sup>***</sup><br>(0.04)</td><td>0.35<sup>***</sup><br>(0.04)</td><td>0.41<sup>***</sup><br>(0.04)</td></tr><tr><td>θ<sub>3</sub></td><td>0.05<sup>***</sup><br>(0.01)</td><td>0.09<sup>***</sup><br>(0.01)</td><td>0.02<br>(0.01)</td></tr><tr><td>θ<sub>4</sub></td><td>-0.06<sup>***</sup><br>(0.01)</td><td>-0.07<sup>***</sup><br>(0.01)</td><td>-0.03<sup>**</sup><br>(0.01)</td></tr><tr><td>θ<sub>5</sub></td><td>0.20<sup>***</sup><br>(0.04)</td><td>0.37<sup>***</sup><br>(0.05)</td><td>0.11<sup>**</sup><br>(0.04)</td></tr><tr><td>θ<sub>6</sub></td><td>0.00<br>(0.00)</td><td>0.00<br>(0.00)</td><td>0.00<br>(0.00)</td></tr><tr><td>θ<sub>7</sub></td><td>-0.00<br>(0.00)</td><td>-0.01<br>(0.00)</td><td>-0.01<br>(0.00)</td></tr><tr><td>θ<sub>8</sub></td><td>-0.66<sup>***</sup><br>(0.12)</td><td>-0.41<sup>***</sup><br>(0.10)</td><td>-0.57<sup>***</sup><br>(0.13)</td></tr><tr><td>θ<sub>9</sub></td><td>-0.14<sup>***</sup><br>(0.03)</td><td>-0.15<sup>***</sup><br>(0.03)</td><td>-0.16<sup>***</sup><br>(0.03)</td></tr><tr><td>θ<sub>10</sub></td><td>0.74<sup>**</sup><br>(0.23)</td><td>-0.24<br>(0.18)</td><td>2.12<sup>***</sup><br>(0.29)</td></tr><tr><td>Adj. R<sup>2</sup></td><td>0.94</td><td>0.93</td><td>0.93</td></tr></tbody></table>

### 4.3.1 Cigarette consumption in the US states

We first consider the case of the US cigarette consumption data. Figure 11 shows the resulting cluster structure obtained with the Lewis et al. (2023) method. Interestingly, we do not find fuzzy units. Moreover, the resulting cluster structure is very different from that obtained with the spatial fuzzy GFE proposed in the paper and shown in Fig. 4. While with the spatial GFE, we obtain a distinction between the west and east of the US, in this case, we get fragmented clusters. For instance, Group 1 includes states such as California and New York, which have relatively high incomes and strong public health policies, but also Alabama and Kentucky, where tobacco farming has historically been significant, and smoking rates are typically higher due to cultural acceptance and less stringent regulation. We highlight that smoking rates tend to be lower in urban, highly educated areas, where health awareness is stronger (Hiscock et al. 2012). But given the partition in Fig. 11, we find that states such as New York, California and Massachusetts that are urbanised and known for high education levels, which correlates with lower smoking rates, are in Group 1 along with states such as Alabama,

<footnote>***p < 0.001; **p < 0.01; *p < 0.05</footnote>
US states: clustering structure

<figure><img src="image_19.png" /><figcaption>**Fig. 11** Clustering structure implied by the fuzzy group fixed-effect estimation without spatial clustering (Lewis et al. 2023)</figcaption></figure>

Kentucky, and Tennessee, which are less urbanised and have lower levels of edu-
cation, which generally lead to higher smoking rates.

In sum, while the partition in Fig. 11 may reveal some patterns, it is not entirely
consistent with known determinants of cigarette consumption, as it groups
together states with significantly different smoking behaviours. An east–west par-
tition would likely be more reasonable, given the distinct economic, cultural and
policy differences that influence cigarette use in these two broad regions of the
US. Therefore, the clustering in Fig. 11 is less explainable compared to that in
Fig. 4.

Table 4 shows the regression results in terms of parameter estimates. Compared to Table 1, in Table 4, we find that the intercepts of the two groups are almost the same and are also equal to the intercept of the full model that does not employ clustering. In Table 1, instead, we find the intercepts to be different, thus suggesting different average cigarette consumption levels over the years in the two clusters, given that the covariates assume zero value. We do not find such a distinction in the model without spatial clustering which thus turns out to be less informative in terms of differences. Moreover, an interesting difference is in the coefficient related to the neighbours spillover θ₄. While in Table 1, we find that this is significant only for the US states in the east, in Table 4, we find it to be significant for both clusters. The effect in Group 1 is larger than in Group 2, where we find most of Midwest states.

Due to the lack of fuzzy units, the partition obtained with the Bonhomme and
Manresa (2015) approach is the same as in Fig. 11. Therefore, also the estimated
parameters are very similar. For this reason, the discussion related to the results with
the Lewis et al. (2023) approach also applies to those obtained with this method.
<table><caption>**Table 4** Estimated parameters with the fuzzy group fixed-effect estimation without spatial clustering (Lewis et al. 2023). Cigarette consumption dataset</caption><thead><tr><th></th><th>Model 1</th><th>Model 2</th><th>Model 3</th></tr></thead><tbody><tr><td>θ<sub>0</sub></td><td>0.17***<br>(0.03)</td><td>0.16***<br>(0.03)</td><td>0.17***<br>(0.03)</td></tr><tr><td>θ<sub>1</sub></td><td>0.96***<br>(0.01)</td><td>0.95***<br>(0.01)</td><td>0.97***<br>(0.01)</td></tr><tr><td>θ<sub>2</sub></td><td>-0.10***<br>(0.01)</td><td>-0.12***<br>(0.01)</td><td>-0.08***<br>(0.01)</td></tr><tr><td>θ<sub>3</sub></td><td>0.03***<br>(0.01)</td><td>0.04***<br>(0.01)</td><td>0.02***<br>(0.01)</td></tr><tr><td>θ<sub>4</sub></td><td>0.04***<br>(0.01)</td><td>0.06***<br>(0.01)</td><td>0.02*<br>(0.01)</td></tr><tr><td>Adj. R<sup>2</sup></td><td>0.97</td><td>0.97</td><td>0.97</td></tr></tbody></table>

### 4.3.2 Insurance demand in the Italian provinces

We then apply the fuzzy GFE approach without spatial dimension to the empirical case of Italian insurance demand. Figure 12 shows the obtained partition, which now does not clearly divide Italy into north and south as the previous results. Mixing provinces from the north with those from the south is not problematic in itself, as some regions in both parts of the country share similar insurance demand characteristics. However, the inconsistencies arise when provinces with vastly different economic and social realities are grouped together, leading to a distorted representation of insurance consumption patterns.

The partition of Italian provinces reveals several inconsistencies, particularly in how provinces with vastly different economic profiles are grouped together. Wealthy, industrialised northern provinces such as Alessandria and Ancona are grouped with economically struggling southern provinces such as Avellino and Agrigento, despite their significant differences in income levels, insurance markets and economic activity. Similarly, economically similar provinces such as Bologna and Modena are placed in different clusters, while Modena is grouped with southern provinces like Catania, which have lower wealth and insurance demand.

Additionally, southern provinces with comparable economic challenges, such as Caserta and Benevento, are split into different groups. Some northern provinces, like Torino, are grouped with southern ones like Reggio Calabria, despite stark differences in economic development and insurance market penetration. Even within northern regions, provinces such as Verona and Vicenza, which share similar industrial profiles, are placed in separate clusters. These inconsistencies result in a partition that does not accurately reflect the determinants of insurance demand across the country.

Table 5 shows the parameter estimates. Several differences can be highlighted between Table 5, showing the results of the GFE method without spatial clustering, and Table 3 showing the results of the proposed approach. In both tables, the

<footnote>***p < 0.001; **p < 0.01; *p < 0.05</footnote>
<figure><img src="image_20.png" /><figcaption>**Fig. 12** Clustering structure implied by the fuzzy group fixed-effect estimation without spatial clustering (Lewis et al. 2023). Insurance consumption dataset</figcaption></figure>

intercept values for all models are negative and significant. However, in Table 3, the differences between the two clusters (north and south) are larger compared to Table 5, suggesting that spatial clustering significantly impacts the baseline insurance demand levels. The coefficients for $\theta_3$, $\theta_4$ and $\theta_5$ differ more noticeably across tables. For instance, $\theta_3$ ranges from 0.05 to 0.09 in Table 3, compared to a narrower range in Table 5. Again, this suggests that spatial clustering introduces more variation in these parameters, reflecting local differences across regions. The greatest difference lies in $\theta_{10}$ about the trust variable. While in Table 3, we found a strong positive coefficient for units in the South cluster, now the coefficient is still significant only for one cluster but its magnitude is much lower. This happens because the partition obtained with the Lewis et al. (2023) approach mixes both northern and southern provinces so that the total effect given by the southern provinces is mitigated by the presence of the northern ones.

Also in this case, the absence of fuzziness allows getting the same results with both Lewis et al. (2023) and Bonhomme and Manresa (2015), and no more differences can be highlighted between the two approaches in terms of partition and results. The results obtained with the proposed spatial GFE allows for highlighting
<table><caption>**Table 5** Estimated parameters with the fuzzy group fixed-effect estimation without spatial clustering (Lewis et al. 2023). Insurance consumption dataset</caption><thead><tr><th></th><th>Model 1</th><th>Model 2</th><th>Model 3</th></tr></thead><tbody><tr><td>θ<sub>0</sub></td><td>-7.54***<br>(0.79)</td><td>-7.96***<br>(0.70)</td><td>-6.78***<br>(0.90)</td></tr><tr><td>θ<sub>1</sub></td><td>0.91***<br>(0.08)</td><td>0.94***<br>(0.08)</td><td>0.86***<br>(0.08)</td></tr><tr><td>θ<sub>2</sub></td><td>0.43***<br>(0.04)</td><td>0.41***<br>(0.05)</td><td>0.48***<br>(0.04)</td></tr><tr><td>θ<sub>3</sub></td><td>0.05***<br>(0.01)</td><td>0.03*<br>(0.01)</td><td>0.07***<br>(0.01)</td></tr><tr><td>θ<sub>4</sub></td><td>-0.06***<br>(0.01)</td><td>-0.04***<br>(0.01)</td><td>-0.07***<br>(0.01)</td></tr><tr><td>θ<sub>5</sub></td><td>0.20***<br>(0.04)</td><td>0.20***<br>(0.04)</td><td>0.17***<br>(0.05)</td></tr><tr><td>θ<sub>6</sub></td><td>0.00<br>(0.00)</td><td>0.00<br>(0.00)</td><td>0.00<br>(0.00)</td></tr><tr><td>θ<sub>7</sub></td><td>-0.00<br>(0.00)</td><td>-0.01*<br>(0.00)</td><td>0.00<br>(0.00)</td></tr><tr><td>θ<sub>8</sub></td><td>-0.66***<br>(0.12)</td><td>-0.60***<br>(0.11)</td><td>-0.90***<br>(0.13)</td></tr><tr><td>θ<sub>9</sub></td><td>-0.14***<br>(0.03)</td><td>-0.16***<br>(0.03)</td><td>-0.13***<br>(0.03)</td></tr><tr><td>θ<sub>10</sub></td><td>0.74**<br>(0.23)</td><td>1.12***<br>(0.23)</td><td>0.20<br>(0.23)</td></tr><tr><td>Adj. R<sup>2</sup></td><td>0.94</td><td>0.93</td><td>0.93</td></tr></tbody></table>

stronger differences within the country in terms of insurance demand compared to non-spatial methods.

# 5 Conclusions

This paper presents an innovative approach to estimating panel data regressions with spatial clustering, highlighting its importance in empirical economics and distinguishing it from previous research.

While conventional cross-sectional regression techniques have proposed methods to improve spatial clustering, the previous studies on panel data have largely neglected this spatial aspect in cluster assignment. Incorporating spatial considerations into clustering enriches our understanding of complex economic dynamics and is consistent with the core tenet of geographical laws. Consequently, this paper makes two important contributions to the existing literature. First, it introduces a refined spatially clustered approach to estimating group-specific panel regression parameters amidst unknown clusters, thereby enhancing precision and robustness.

<footnote>***p < 0.001; **p < 0.01; *p < 0.05</footnote>
Second, it pioneers the use of fuzzy logic in the estimation of spatially clustered panel data models.

To validate the efficacy of the proposed method, two empirical case studies are presented, drawing on different spatial datasets from the existing literature. In the first case, examining cigarette consumption across US states (Baltagi and Levin 1992), we identify statistically significant effects of neighbouring state cigarette prices only within a cluster of western US states, in contrast with the broader US findings in Baltagi and Levin (1992). In the second case, analysing non-life insurance consumption in Italian provinces (Millo and Carmeci 2011), we uncover differences in determinants and their magnitudes between northern and southern Italy. In particular, variables such as trust in others, significant according to Millo and Carmeci (2011), show statistically significant effects in the south but not in the north, demonstrating regional differences in economic dynamics.

The fact the clusters align with intuitive geographical or economic divisions of both the US and Italy is a strength of the proposed model. Indeed, this alignment shows that the clustering results obtained with the spatial GFE are not arbitrary but are driven by meaningful and well-established economic and spatial factors. Since the spatial GFE clustering approach naturally reconstructs these known spatial patterns, this validates the model’s ability to detect meaningful distinctions. However, while general knowledge may suggest certain divisions, such as those based on geography, we highlight that the spatial GFE model goes beyond this by not only confirming these divisions but also providing a data-driven explanation for such. This ensures that the clusters are not simply the result of preconceived assumptions, but are grounded in the underlying data structure.

Future research should focus on developing models that properly address endogeneity, a still overlooked aspect of clustered regression models that are needed to improve the reliability of findings from a causal perspective. Techniques based on clustered instrumental variable estimation, for example, could help to mitigate this problem. However, dealing with endogeneity in spatial contexts poses unique challenges (e.g. see Betz et al. 2018) that need to be taken into account. Future research should focus on innovative methods to disentangle spatial effects from endogeneity in contexts where unknown clustered effects exist.

## References

Accetturo, A., Albanese, G., De Paola, M., et al.: The north-south gap: economic development and public intervention. Ital. Econ. J. 10, 957-975 (2024)

Ahearn, B.: Finance-led divergence in the regions of Italy. Financ. Hist. Rev. 12(1), 7-41 (2005)

Alciati, M.H., Frosh, M., Green, S.B., et al.: State laws on youth access to tobacco in the united states: measuring their extensiveness with a new rating system. Tob. Control 7(4), 345-352 (1998)

Alin, A.: Simpson's paradox. Wiley Interdiscip. Rev.: Comput. Stat. 2(2), 247-250 (2010)

Anselin, L., Arribas-Bel, D.: Spatial fixed effects and spatial dependence in a single cross-section. Pap. Reg. Sci. 92(1), 3-17 (2013)

Aquaro, M., Bailey, N., Pesaran, M.H.: Estimation and inference for spatial models with heterogeneous coefficients: an application to us house prices. J. Appl. Econom. 36(1), 18-44 (2021)

Baltagi, B.H., Levin, D.: Estimating dynamic demand for cigarettes using panel data: the effects of bootlegging, taxation and advertising reconsidered. Rev. Econ. Stat. pp 148-155 (1986)
Baltagi, B.H., Levin, D.: Cigarette taxation: raising revenues and reducing consumption. *Struct. Chang. Econ. Dyn.* **3**(2), 321–335 (1992)

Bester, C.A., Hansen, C.B.: Grouped effects estimators in fixed effects models. *J. Econom.* **190**(1), 197–208 (2016)

Betz, T., Cook, S.J., Hollenbach, F.M.: On the use and abuse of spatial instruments. *Polit. Anal.* **26**(4), 474–479 (2018)

Bezdek, J.C., Ehrlich, R., Full, W.: FCM: the fuzzy c-means clustering algorithm. *Comput. Geosci.* **10**(2–3), 191–203 (1984)

Billé, A.G., Benedetti, R., Postiglione, P.: A two-step approach to account for unobserved spatial heterogeneity. *Spat. Econ. Anal.* **12**(4), 452–471 (2017)

Bonhomme, S., Manresa, E.: Grouped patterns of heterogeneity in panel data. *Econometrica* **83**(3), 1147–1184 (2015)

Botello-Harbaum, M.T., Haynie, D.L., Iannotti, R.J., et al.: Tobacco control policy and adolescent cigarette smoking status in the united states. *Nicotine Tob. Res.* **11**(7), 875–885 (2009)

Brusco, M.J., Cradit, J.D., Steinley, D., et al.: Cautionary remarks on the use of clusterwise regression. *Multivar. Behav. Res.* **43**(1), 29–49 (2008)

Coppi, R., D'Urso, P., Giordani, P., et al.: Least squares estimation of a linear regression model with LR fuzzy response. *Comput. Stat. Data Anal.* **51**(1), 267–286 (2006)

Coppi, R., D'Urso, P., Giordani, P.: A fuzzy clustering model for multivariate spatial time series. *J. Classif.* **27**, 54–88 (2010)

Dembele, D., Kastner, P.: Fuzzy c-means method for clustering microarray data. *Bioinformatics* **19**(8), 973–980 (2003)

DeSarbo, W.S., Cron, W.L.: A maximum likelihood methodology for clusterwise linear regression. *J. Classif.* **5**, 249–282 (1988)

Di Mari, R., Rocci, R., Gattone, S.A.: Lasso-penalized clusterwise linear regression modelling: a two-step approach. *J. Stat. Comput. Simul.* **93**(18), 3235–3258 (2023)

D'Urso, P., Santoro, A.: Fuzzy clusterwise linear regression analysis with symmetrical fuzzy output variable. *Comput. Stat. Data Anal.* **51**(1), 287–313 (2006)

D'Urso, P., Massari, R., Santoro, A.: A class of fuzzy clusterwise regression models. *Inf. Sci.* **180**(24), 4737–4762 (2010)

D'Urso, P., De Giovanni, L., Disegna, M., et al.: Fuzzy clustering with spatial-temporal information. *Stat. Stat.* **30**, 71–102 (2019)

D'Urso, P., De Giovanni, L., Vitale, V.: Spatial robust fuzzy clustering of COVID 19 time series based on b-splines. *Spat. Stat.* **49**(100), 518 (2022)

D'Urso, P., De Giovanni, L., Alaimo, L.S., et al.: Fuzzy clustering with entropy regularization for interval-valued data with an application to scientific journal citations. *Ann. Oper. Res.* **342**, 1605–1628 (2023)

D'Urso, P., et al.: Fuzzy clustering. *Handbook of cluster analysis*, Chapman and Hall/CRC, Boca Raton, pp 545–574 (2015)

Elhorst, J.P.: Specification and estimation of spatial panel data models. *Int. Reg. Sci. Rev.* **26**(3), 244–268 (2003)

Fouedjio, F.: A hierarchical clustering method for multivariate geostatistical data. *Spat. Stat.* **18**, 333–351 (2016)

Gruber, J., Poterba, J.: Tax incentives and the decision to purchase health insurance: evidence from the self-employed. *Q. J. Econ.* **109**(3), 701–733 (1994)

Guiso, L., Paiella, M.: Risk aversion, wealth, and background risk. *J. Eur. Econ. Assoc.* **6**(6), 1109–1150 (2008)

Hathaway, R.J., Bezdek, J.C.: Switching regression models and fuzzy clustering. *IEEE Trans. Fuzzy Syst.* **1**(3), 195–204 (1993)

Hennig, C.: Identifiability of models for clusterwise linear regression. *J. Classif.* **17**, 273–296 (2000)

Hiscock, R., Bauld, L., Amos, A., et al.: Socioeconomic status and smoking: a review. *Ann. N. Y. Acad. Sci.* **1248**(1), 107–123 (2012)

Lawrence, D., Fagan, P., Backinger, C.L., et al.: Cigarette smoking patterns among young adults aged 18–24 years in the united states. *Nicot. Tob. Res.* **9**(6), 687–697 (2007)

Lee, L.-f., Yu, J.: Estimation of spatial autoregressive panel data models with fixed effects. *J. Econom.* **154**(2), 165–185 (2010)

Lee, J., Gangnon, R.E., Zhu, J.: Cluster detection of spatial regression coefficients. *Stat. Med.* **36**(7), 1118–1133 (2017)
Lewis, D.J., Melcangi, D., Pilossoph, L., et al.: Approximating grouped fixed effects estimation via fuzzy clustering regression. J. Appl. Econom. 38(7), 1077–1084 (2023)

Li, D., Moshirian, F., Nguyen, P., et al.: The demand for life insurance in OECD countries. J. Risk Insur. 74(3), 637–652 (2007)

Mattera, R.: A weighted approach for spatio-temporal clustering of COVID-19 spread in Italy. Spat. Spa-tio-tempor. Epidemiol. 41(100), 500 (2022)

Mattera, R., Franses, P.H.: Are African business cycles synchronized? evidence from spatio-temporal modeling. Econ. Model. 128(106), 485 (2023)

Menon, M., Perali, F., Ray, R., et al.: Heterogeneity in prices and cost of living within a country: new evidence on the north-south divide in Italy. Econ. Model. 126, 106401 (2023)

Millo, G.: The income elasticity of nonlife insurance: a reassessment. J. Risk Insur. 83(2), 335–362 (2016)

Millo, G., Carmeci, G.: Non-life insurance consumption in Italy: a sub-regional panel data analysis. J. Geogr. Syst. 13, 273–298 (2011)

Nelson, G.D., Rae, A.: An economic geography of the United States: from commutes to megaregions. PLoS ONE 11(11), e0166,083 (2016)

Outreville, J.F.: The economic significance of insurance markets in developing countries. J. Risk Insur. 57(3), 487–498 (1990)

Pal, N.R., Bezdek, J.C.: On cluster validity for the fuzzy c-means model. IEEE Trans. Fuzzy Syst. 3(3), 370–379 (1995)

Pan, T., Hu, G., Shen, W.: Identifying latent groups in spatial panel data using a Markov random field constrained product partition model (2020). arXiv preprint arXiv:2012.10541

Pham, D.L.: Spatial models for fuzzy clustering. Comput. Vis. Image Underst. 84(2), 285–297 (2001)

Potts, R.B.: Some generalized order-disorder transformations. In: Mathematical proceedings of the cam-bridge philosophical society. Cambridge University Press, Cambridge, pp 106–109 (1952)

Salvati, L.: ‘North of the south’ or ‘south of the north’? Revisiting the spatially complex economic divide in Italy. J. Urban and Regional Anal, 5(1), 65–75 (2013)

Späth, H.: Algorithmus 39. klassenweise lineare regression. Computing 22, 367–373 (1979)

Sugasawa, S.: Grouped heterogeneous mixture modeling for clustered data. J. Am. Stat. Assoc. 116(534), 999–1010 (2021)

Sugasawa, S., Murakami, D.: Spatially clustered regression. Spat. Stat. 44(100), 525 (2021)

Tobler, W.R.: A computer movie simulating urban growth in the detroit region. Econ. Geogr. 46(sup1), 234–240 (1970)

Von Lyncker, K., Thoennessen, R.: Regional club convergence in the EU: evidence from a panel data analysis. Empir. Econ. 52, 525–553 (2017)

Wang, X.Y., Garibaldi, J.M., Bird, B., et al.: Novel developments in fuzzy clustering for the classification of cancerous cells using FTIR spectroscopy. In: Pedrycz, W. (ed.) Advances in fuzzy clustering and its applications, p. 405. Wiley, New Jersey (2007)

Yu, J., De Jong, R., Lf, Lee: Quasi-maximum likelihood estimators for spatial dynamic panel data with fixed effects when both n and t are large. J. Econom. 146(1), 118–134 (2008)

Yu, D., Zhang, Y., Wu, X., et al.: The varying effects of accessing high-speed rail system on China's county development: a geographically weighted panel regression analysis. Land Use Policy 100(104), 935 (2021)

Zhu, H., Fan, J., Kong, L.: Spatially varying coefficient model for neuroimaging data with jump discontinuities. J. Am. Stat. Assoc. 109(507), 1084–1098 (2014)

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.
## Authors and Affiliations

Roy Cequeti¹,² · Pierpaolo D'Urso¹ · Raffaele Mattera¹,³

<figure><img src="image_21.png" /></figure>

Raffaele Mattera
raffaele.mattera@uniroma1.it

Roy Cequeti
roy.cerqueti@uniroma1.it

Pierpaolo D'Urso
pierpaolo.durso@uniroma1.it

¹ Department of Social and Economic Sciences, Sapienza University of Rome, Piazzale Aldo Moro, 5, 00185 Rome, Italy

² GRANEM, University of Angers, Angers, France

³ Department of Mathematics and Physics, University of Campania “Luigi Vanvitelli”, Viale Abramo Lincoln, 5, 81100 Caserta, Italy