Obesity and labour market outcomes in Italy: a dynamic panel data evidence with correlated random effects.

This paper investigates the effects of obesity, socio-economic variables, and individual-specific factors on work productivity across Italian regions. A dynamic panel data with correlated random effects is used to jointly deal with incidental parameters, endogeneity issues, and functional forms of misspecification. Methodologically, a hierarchical semiparametric Bayesian approach is involved in shrinking high dimensional model classes, and then obtaining a subset of potential predictors affecting outcomes. Monte Carlo designs are addressed to construct exact posterior distributions and then perform accurate forecasts. Cross-sectional Heterogeneity is modelled nonparametrically allowing for correlation between heterogeneous parameters and initial conditions as well as individual-specific regressors. Prevention policies and strategies to handle health and labour market prospects are also discussed.

Introduction

This paper aims to improve the existing literature on the relationships between obesity, healthcare, and labour market outcomes using a dynamic panel data (DPD) with correlated random effects (CRE). The main thrust of this study is to examine in depth these relationships across Italian regions using: (i) additional time-varying socio-economic indicators and exogenous factors to deal with endogeneity issues (because of unobserved heterogeneity and omitted variables) and functional forms of misspecification; (ii) unobserved time-invariant individual-specific heterogeneity affecting outcomes; (iii) semiparametric Bayesian inference to avoid the curse of dimensionality and variable selection problems in high dimensional settings; and (iv) conditional density forecasts performing policy-relevant strategies.

Obesity is a complex condition increasing the risk of various health problems, such as cancer, diabetes, hypertension, and depression, which may negatively affect economic outcomes and thus the probability of employment (see, e.g., [13, 23, 58, 67], and [4]). Furthermore, employers may discriminate against obese workers, resulting in worse labor market performance (see, e.g., [11, 36], and [51]). In that context, it is not always clear as discrimination may be directly related to labour market outcomes such as social participation, per-hour wage rates, sales (see, e.g., [17] and [9]). Earlier studies highlight that physical appearance may be directly related to productivity, affecting the individual’s capacity to work. For instance, obese people may be related to non-desirable personality traits potentially affecting productivity, becoming less socially and intellectually skilled than their non-obese counterparts (see, e.g., [60, 61, 66], and [69]).

The previous literature has primarily looked at the effects of obesity on wages that would be not necessarily the same as the ones of obesity on employment. However, if wages are not fully flexible or obesity causes serious health problems, there would be a severe causal-effect of obesity on employment (see, for instance, [1]). In this context, investigating the effect of obesity on labour market outcomes is complicated due to potential reversed causality and endogeneity issues. Moreover, there may also be unobserved heterogeneous characteristics that vary systematically between obese and non-obese people, affecting employment. For instance, whether high discount rates can encourage people to be more prone to weight-gaining consumption, they would make investments in human capital and thus future labour market outcomes less attractive ( [20, 21] and [22]).

Important issues of the effects of obesity on labour market outcomes can be further addressed focusing on groups having traditionally lower participation rates such as females. The existing literature on the effects between obesity and wage productivity is still limited, reaching mixed conclusions. In [20, 21] and [55], no effect of obesity on employment among US men and women is obtained. On the other hand, in [53, 54], a significant and negative effect for both genders in the UK are found. Different results are found in Italy instead. Even if the country is characterized by one of the slimmest adult populations in Europe (OECD 20171), Italy is experiencing a strong increase in the portion of obese people, mainly among children age, making it a particularly relevant problem to be addressed (see, e.g., [16]). Moreover, because Italy has always been the country of fashion, larger attention is devoted to physical appearance. Thus, according to international studies, weight-based discrimination would be a relevant topic, strongly affecting every aspect of obese people’s everyday life including labour market outcomes (see, e.g., [18] and [65]). However, since obesity denotes a major issue in every developed country, the existing literature has devoted much less attention to the weight bias in the job market compared to race and gender discrimination (see, for instance, [2]).

The relationship between high body weight and labour market outcomes has been primarily studied using data from developed and high income countries, such as the US and West Europe (e.g., England, Denmark, and Finland). In the last decade, the main studies concerning labour market outcomes have focused on wages/earnings, employment, and occupational selection (see, e.g., [28, 14, 15, 44, 62, 29, 42], and [30]). Related papers, focused on the US, have used the National Longitudinal Survey of Youth (NLSY) data, founding mixed results (see, e.g., [63, 50], and [57]). Some shortcomings of these studies are that they ignore potential endogeneity and spatial interdependence of obesity accounting for small and unrepresentative samples and cross-sectional data (see, for instance, [26] and [34]).

Concerning DPD models, they are widely used in empirical economics for forecasting individuals’ future outcomes (see, e.g., [33, 38, 46], and [48]), and allowing the possibility of controlling for unobserved time-invariant individual heterogeneity (see, e.g., [11] and [7] for linear case). Such heterogeneity is an important issue and failure to control for it results in misleading inferences. That problem is even more severe when the unobserved heterogeneity may be correlated with covariates. However, in these traditional methods, with large cross-sectional dimension (N) and short (fixed) time-series (T), the estimators of homogenous parameters would result biased and inconsistent due to incidental parameter problems (because of individual-specific effects). Indeed, leaving the individual heterogeneity unrestricted, the number of individual-specific effects would grow with the sample size. Moreover, with short T, the estimates of the heterogeneous parameters would be highly contaminated from idiosyncratic shocks obtaining inaccurate forecasts.

The methodology proposed in this paper focuses on the aforementioned issues and consists of a hierarchical semiparametric Bayesian approach for inference in dynamic panel setup with cross-sectional heterogeneity. Methodologically, Conjugate Informative Proper Mixture (CIPM) priors are used to obtain a subset of best combination of predictors (or model solutions) affecting outcomes of interest, where best stands for the model solution providing the most accurate predictive performance over all candidate models. Markov Chain Monte Carlo (MCMC) algorithms and implementations are used to construct exact posterior distributions, and then perform better density forecasts and policy issues.

The contributions of this paper are threefold. First, full Bayesian methods are used to deal with potential features in real-world data such as density forecasts, incidental parameters, endogeneity issues, and structural model uncertainty (because of one or more parameters are posited as the source of model misspecification problems). CIPM priors and MCMC-based Posterior Model Probabilities (PMPs) are used to: (i) include cross-sectional individual information from the whole panel, acting as a strong model selection in high dimensional model classes; (ii) outline assumptions and specification strategy of the informative priors on covariates, outcomes of interest, and the distribution of unobserved heterogeneity; and (iii) deal with overfitting and model uncertainty when addressing variable selection problems. Here, ’overfitting’ refers to the overestimation of effect sizes since more complex models (or models with a larger number of unknown parameters) always provide a somewhat better fit to the data than simpler models.

Second, I build on and extend the [56]’s analysis who develops a Robust Open Bayesian (ROB) procedure in two stages for implementing Bayesian Model Averaging (BMA) and Bayesian Model Selection (BMS) in multiple high dimensional linear regression models. More precisely, I extend the prior specification strategy in a panel context with either dynamic or time-invariant factors, taking the name of Time-varying ROB (TROB) procedure. The main novelty of the implementation is the use of nonparametric Bayesian inference for heterogeneous parameters. Conversely, homogeneous parameters are modelled through parametric priors. In the first stage, a subset of covariates better explaining and thus fitting the data are chosen according to their PMPs. The latter denotes the probability of each candidate model performing the data. In the second stage, further shrinkage is performed to obtain a smallest final subset of best submodels having statistically significant predictive capability. Finally, the submodel with the highest Bayes Factor in logarithm (lBF) would correspond to the final solution containing a potential subset of candidate predictors with a higher significant overall F value and sufficiently strong adjusted- () measure, where ’strong’ refers to value equal to or bigger than .

Third, MCMC algorithms and implementations are addressed to construct exact posterior distributions dealing with semiparametric forecasting problems. Better density forecasts are involved in DPB with CREs because of three main features: (i) the use of a hierarchical full Bayesian approach; (ii) the framework of informative mixture priors; and (iii) the observation of incidental parameters.

An empirical example across Italian regions describes the estimating procedure and forecasting performance. It covers the period , where conditional density forecasts are performed extending it up to 2021 to also address potential findings concerning the impact of the ongoing pandemic disease on the global economy. Then, the sample is divided into three groups to better highlight the strengths and limits of policy strategies: (1) Northern Italy; (2) Central Italy; and (3) Southern Italy. In this study, natural conjugate priors are used to deal with persistent structural breaks affecting the panel data analysis: 2007, due to the global financial crisis; 2011, due to fiscal consolidations in post-financial crisis; and 2014, due to recovery fund.

The remainder of this paper is organized as follows. Section "Econometric method and estimation procedure" introduces the econometric model, the Bayesian inference, and the estimating procedure. Section "Multivariate dynamic panel data" describes the dynamic procedure to outline prior specification strategy, and then construct posterior distributions. Section "Bayesian analysis" displays theoretical properties concerning dynamic panel data and the semiparametric Bayesian approach. Section "Dynamic analysis" describes the data, the empirical analysis, and density forecasts along with policy issues. The final section contains some concluding remarks.

Econometric method and estimation procedure

Multivariate dynamic panel data

The baseline dynamic panel data (DPD) model with correlated random effects (CREs) and cross-sectional heterogeneity is:where i = 1, 2, ..., N, t = 1, 2, ..., T, is a vector of outcomes, and are vectors of predetermined and endogenous variables for each i, respectively, with and denoting autoregressive coefficients to be estimated for each i, is a vector of strictly exogenous factors for each i, with denoting the regression coefficients to be estimated, is a heterogeneous intercept containing, for example, time-constant differences among regions, and is a vector of unpredictable shock (or idiosyncratic error term), with and if and , and otherwise.

In this study, for notational simplicity, I assume only one lag for every time-varying parameters. Here, some considerations are in order: (i) the predetermined variables refer to the lagged outcomes () containing, for example, control variables; (ii) the ’s denote cross-sectional heterogeneity affecting the outcomes ; (iii) CREs matter and then the ’s are treated as random variables and possibly correlated with some of the covariates within the panel; (iv) the strictly exogenous factors contain dummy variables to test, for example, the presence of policy shifts (or structural breaks); and (v) the instruments are fitted values from AutoRegressive (AR) parameters based on time-varying lagged variables.

Furthermore, because of the location aspects in Eq. (1), two potential problems estimating the regression coefficients in Eq. (1) need to be accounted for: spatial dependence among observations in the system and spatial heterogeneity among time-varying parameters in the relationships. Both of them lead to violations of assumptions underlying unbiasedness, consistency, and efficiency (see, for instance, [10] and [40]). According to Eq. (1), spatial interdependence is involved in the estimation including AR coefficients, and spatial heterogeneity is addressed let the model be multidimensional (panel data).

Assuming stationarity to hold in Eq. (1), the time-series regressions are valid and the estimates feasible. Thus, some moment restrictions need to hold to address exact identification in a context of CREs, and then estimate the time-varying parameters ( and ) for (see, for instance, [3, 8], and [12]). More precisely, I assume that and are independently distributed across i and have the familiar error components structure:andThen, I also assume the standard assumption concerning the initial conditions :

Bayesian analysis

The main thrust of the TROB procedure is threefold. First, it performs a strong model selection to obtain the best model solution better explaining and thus fitting the data in high dimensional panel setups. Second, CIPM priors allow for modelling nonparametrically common (, , ) and heterogeneous (, ) coefficients, and then their distribution changing in a corresponding fashion in accordance with different model solutions. Third, MCMC algorithms and implementations can be constructed to perform better density forecasts concerning (potential) prevention policies and strategies when investigating health and labour market prospects.

Given the DPD in Eq. (1), I decompose the vectors of the directly observed (endogenous) variables: , with denoting lagged outcomes to capture the persistence and including lagged control variables such as individual-specific factors; and , with referring to economic conditions and standing for health and social prevention activities. Then, I combine the homogeneous parameters into the vector .

To model the key latent heterogeneities () and observed determinants (, , ) when dealing with high dimensional settings, I define the conditioning set at period t () and the structural density () as:andThe error terms () are individual-time-specific shocks characterized by zero mean and homoskedastic Gaussian innovations. In a hierarchical framework, I combine the individual-specific heterogeneity into the vector under cross-sectional heterogeneity and homoskedasticity. Assuming correlated random coefficients model, and may be correlated with each other, with:Given these specifications, the DPD model in Eq. (1) tends to be less parsimonious and suitable for density forecasting. Thus, I perform a fully Bayesian approach to estimate the prior distribution of the correlated random effects by pooling the cross-sectional information from the whole panel. The main difference between an empirical and fully Bayesian approach is that the former picks the ’s distribution by maximizing the maximum likelihood of the data (see, e.g., [25, 38, 45, 41], and [32, 33]). Conversely, fully Bayesian methods construct a prior for the correlated random effects and then evaluates it in view of the observed panel data (see, for instance, [46] and [48] (linear case); and [47] (non-linear case)). In addition, the fully Bayesian approach tends to be more suitable for density forecasting and more easily for extending the TROB procedure to the nonlinear case.

Let be the full panel set containing all (potential) model solutions, the variable selection problem is addressed by imposing an auxiliary indicator variable , with , containing every possible subset choices, where if is small (absence of j-th covariate in the model) and if is sufficiently large (presence of j-th covariate in the model). According to the [56]’s framework, I run the TROB procedure by pooling all high dimensional model solutions () to shrink the model and the parameter space at the same time. The shrinking process jointly deals with overfitting (or overestimation of effect sizes) and model uncertainty (involved in the procedure) by Posterior Model Probabilities (PMPs) for every candidate models. They are defined as:where denotes the multidimensional (natural) parameter space for , denotes a countable collection of all (potential) model solutions for any individual i. The integrand in Eq. (8) is defined as:where denotes the joint likelihood and is the marginal likelihood, with referring to the conditional prior distribution of and .

Let N be high dimensional and T be sufficiently large, the calculation of the integral is unfeasible, and then MCMC integrations are needed to move through the model and the parameter space at the same time to obtain a reduced set containing the best combination of predictors. Thus, I perform the two steps involved in the TROB procedure dealing with high multidimensional dynamic settings. The first stage entails obtaining a subset which contains the best model solutions:where denotes the submodel solutions of the DPD in Eq. (1), with , , , and is a threshold chosen arbitrarily for an enough posterior consistency2. Because of the panel considered in this study features high dimensional cross-sectional units N (predictors ) and non-fixed time-series T, I use to jointly manage all equations within the system (through the conditioning set ), and their potential interactions (through time-varying common covariates).

The second stage entails reducing the model space and then obtaining a smaller subset of best submodel solutions:where , denotes the PMPs, with referring to a new auxiliary variable containing the only best model solutions in the subset and referring to a new arbitrary threshold to evaluate the probability of the model solutions in performing the data (PMPs). In this study, I still use for , independently of N, for a sufficient prediction accuracy in explaining the data.

The final model solution () consists of a subset of the submodel solutions obtained in the second stage () to perform conditional density forecasting and policy-making. It will correspond to the outcome with higher log natural Bayes Factor (lBF):In this study, the lBF is interpreted according to the scale evidence in [56], but with more stringent conditions:

Dynamic analysis

Prior specification strategy and Bayesian predictor

The variable selection procedure entails estimating and as posterior means (the probability that a variable is in the model). All observed variables in and individual heterogeneity in are hierarchically modelled via CIPM priors:wherewith and standing for Normal and independent Inverse-Gamma distributions, respectively, referring to the cross-sectional information available at time , and denoting the model prior choice related to the sum of the PMPs (or Prior Inclusion Probabilities) with respect to the model size , through which the ’s will require a non-0 estimate or the ’s should be included in the model. In this way, one would weigh more according to model size and, setting large for smaller , assign more weight to parsimonious models. Moreover, the use of independent distribution allows cross-equation independence of the coefficient distributions and remove the dependence of and (see, for instance, [27]).

All hyperparameters are known. More precisely, collecting them in a vector , where , they are treated as fixed and are either obtained from the data to tune the prior to the specific applications (such as ) or selected a priori to produce relatively loose priors (such as ). Here, is restricted to a benchmark prior according to the non-0 components of .

Nevertheless, to accomodate the correlated random coefficients model where the individual-specific heterogeneity () can be correlated with the conditioning variables and then , I use a nonparametric Bayesian prior accomodating conditional density forecasting. As suggested in [59], I consider the Mixture of Gaussian Linear Regressions (MGLR) prior and extend it to the multidimensional setup conditioning on :where is a component label, are component parameters, with denoting a conjugate multivariate-normal-inverse-gamma (or multivariate-normal-inverse-Wishart) distribution for scalar , and is a component probability. The mixture probabilities are characterized by a Probit Stick-Breaking Process:where is a common prior, with as in [64]. This prior allows to borrow information across populations by shrinking the stick-breaking ratios that corresponds to the l-th mixture component towards a common value, with l denoting the distribution of . However, this formulation would hold for random effect models, and then I implement it by evaluating as a stochastic function drawn from Gaussian process with zero mean and variance , where in homoskedastic case. In this way, I am able to estimate it by drawing posterior distributions for .

With these specifications, four considerations are in order: (i) component means are linear in ; (ii) component covariances are independent of ; (iii) mixture probabilities can be modelled as functions of ; and (iv) CIPM priors involved in TROB procedure absorb dependency on , without requiring any further dependency of component means and covariances on beyond the MGLR specification. This latter consists of a methodological contribution with respect to the [59]’s framework according to high dimensional panel settings.

In the compound decision theory, the infeasible optimal oracle (or benchmark) predictor implies that and the distribution of the unobserved heterogeneity () are known and the values of are unknown across units. Thus, I construct a feasible semiparametric Bayesian predictor to provide an estimate of this distribution and then perform better conditional density forecasts asymptotically converging to the oracle forecast. The density forecasts of are defined as:where denotes the predictor conditional on the common () and heterogeneous () parameters, with referring to the individual i’s data. In Eq. (21), the first term captures individual i’s uncertainty due to future shock , with , and the second term accounts for individual-specific posteriors corresponding to:with . By construction, dealing with high dimensional time-series, I am able to extend the density forecasts in Eq. (21) to -period-ahead forecasts by pooling the information from the whole panel. In this way, I am able to update all the prior beliefs to construct posterior distributions of , and then the feasible semiparametric Bayesian predictor for -period-ahead forecasts by integrating the conditional predictor over the posteriors. This latter is defined as:where, according to the main goal of this study, the first term reflects functional forms of misspecification due to future shocks on multiple-period-ahead forecasts, and the second term deals with endogeneity issues due to cross-sectional individual-specific heterogeneity in high dimensional dynamic panel setups. Methodologically, the infeasible oracle predictor (or benchmark forecast) would correspond to the semiparametric Bayesian predictor in Eq. (23) by plugging the true (expected) values of the common parameters ) and underlying distribution of () into the conditional predictor in Eq. (21):

MCMC algorithms and density forecasts

Despite the dramatic parameter reduction involved in the TROB procedure, the analytical computation of posterior distributions of is unfeasible, where denotes the expectations of outcomes associated with density forecasts [Eq. (21)]. Thus, I use MCMC implemetations to draw conditional posterior distributions of . More precisely, I include a variant of the Gibbs sampler approach, the Kalman-Filter technique, to analytically obtain forward recursions for posterior means and covariances on . They are, respectively:withThus, the marginal distributions of can be computed by averaging over draws in the nuisance dimensions, and the Kalman filter backward can be run to obtain their posterior distributions:withwhere and are smoothed -period-ahead forecasts of and covariances of the forecast error. The above output of the Kalman filter is used to generate a random trajectory for by the backward recursions on , starting with a draw of from 3.

Given Eq. (28), the other posterior distributions can be defined as:Here, some considerations are in order. In Eq. (31), , with and denoting the arbitrary scale parameter and the posterior predictive variance of according to its distribution (l) in homoskedastic case, respectively. In this analysis, and according to the sample size as described in Eq. (13).

In Eq. (32), refers to the model posterior choice according to the sum of the PMPs obtained in the second stage with respect to model size , with accounting for the non-0 components in .

In Eq. (33), and , with and denoting the arbitrary degrees of freedom (sufficiently small) and the arbitrary scale parameter, respectively. In this analysis, and .

Finally, the last two common parameters to be defined are and . In conjugate settings and flexible mixture probabilities involved in CIPM priors, I use another MCMC integration, the Gibbs sampling, to accommodate the underlying distribution (l) conditional on . Methodologically, this framework would be similar to approximating the conditional density of individual-specific heterogeneity through Bayes’ theorem, but without explicitly modelling the distribution of . In this way, the resulting feasible predictor in (23) would be still close to the oracle one, and then both the estimated common parameters () and the estimated distribution of the cross-sectional heterogeneous effects () achieve posterior consistency. According to the MGLR prior in Eq. (19), the posteriors are:where , , and , with denoting the arbitrary degrees of freedom (sufficiently small), and and referring to the arbitrary scale parameters. In this analysis, , , and .

The accuracy of the density forecasts is measured by the Predictive Score in logarithm (lPS) as in [31]:where denotes the predictive likelihood with respect to the estimated model conditional on the observed data .

Theoretical properties

The proof of the posterior consistency is as follow.

Assumption 4.1

(Identification: General Model) Consider the DPD model in Eq. (1):

Model Setup

() are i.i.d. across i.

For all t, conditional on (), is independent of .

is independent of .

Conditioning on () and , they are independent of each other.

Let is i.i.d. across i and independent of ().

Identification

The characteristic functions for and do not steadily disappear altogether into high-shrinkage processing method.

Let be the composite error at time t, the sequence is almost certainly serially correlated, and definitely is if is serially uncorrelated.

With the panel setup in (1)m, with large N and sufficiently large T, includes interactions of variables with time periods dummies and general non-linear functions and interactions, so the model is quite flexible.

For the CRE approach, each kind of covariates in is separated out and the heterogeneous factors in are correlated with them.

Assumption 4.2

(Identification: (Un)balanced Panel) For all i:

Identification

is observed.

() are observed.

Sequential Exogeneity (Conditional on the Unobserved Effects)

for .

Sequential exogeneity is a middle ground between contemporaneous and strict exogeneity. It allows lagged dependent variables and other variables that change in reaction to past shocks.

Because including contemporaneous exogeneity, standard kinds of endogeneity, where some elements of are correlated with , are ruled out such as measurement error, simultaneity, and time-varying omitted variables.

Sequential exogeneity is less restrictive than strict exogeneity imposing restrictions on economic conditions.

Model Setup

The term “correlated random effects” is used to denote situations where one models the relationship between and .

The CRE approach allows to unify fixed and random effects estimation approaches.

With the CRE approach, time-constant variables can be included within the system.

Statement 4.2.1 (Posterior Distributions) Under Assumptions (4.1) and (4.2), all posterior distributions in (28)–(33) hold and are estimable through MCMC algorithms and implementations.

Theorem 4.3

(Posterior Consistency: Correlated Random Coefficients) Given the DPD in (1):

Model: Assumptions (4.1) and (4.2).

Covariates: satisfy Assumptions (4.1) and (4.2).

Common Parameters:

is unknown and estimable.

The domain of is finite and estimable (homoskedastic case).

is in the interior of supp().

For homoskedastic models, the domain of is bounded by for some , where () is some small (large) positive number.

Stick-breaking Process: The base distribution of the MGLR prior in (19) is characterized by a multivariate normal distribution on and an inverse Wishart distribution on with arbitrary scale parameters sufficiently small and arbitrary degrees of freedom close to zero. Thus, posterior distributions are strongly consistent at .

Statement 4.3.1 (Density Forecasts: Correlated Random Coefficients) Given conditions in Theorem (4.3), density forecasts converge to the oracle.

Theorem 4.4

(Density Forecasts: General Semiparametric Model) Given the DPD in Eq. (1):

Posterior Consistency: Conditions in Theorem (4.3).

Model:

is bounded above by some positive numbers L, with .

The second moment of exists and is positive.

If l is a conditional distribution, there exists such that for all , with denoting an induced measure.

Likelihood:

is continuous in .

, for some .

There exists such that for all and for some :

.

.

Then, given individuals i’s, density forecasts converge to the oracle in probability for all and :

Empirical evidence

Data description and results

The dataset contains 53 observable variables accounting for (potential) covariates described through the vectors , , and . In this study, I split them in six groups: (i) lagged outcomes (); (ii) individual-specific characteristics (); (iii) economic indicators (); (iv) healthcare and social factors (); (v) exogenous variables (); and (vi) lagged predictors concerning some spatial aspects based on the arbitrary delineation of the spatial units (e.g., county boundaries), the spatial distribution of health systems (e.g., health care expenditures), and the inherent mobility of county residents (e.g., physical activity). In this study, spatial units are dealth with the panel framework involved in the estimation, the added value per employee denotes the variable of interest accounted for evaluating cross-unit wage effects, and individual’s weight is used to estimate cross-unit excess weight across Italian regions. The latter has been computed according to the CDC4 as , where kg is a person’s weight in kilograms and is its height in metres squared. It is so defined: underweight (below 18.5); normal (18.5 - 24.9); overweight (25.0 - 29.9); and obese (30.0 and above). In this study, I focus on the last two indicators.

The estimation sample is expressed in years and accounts for all the 21 Italian regions spanning the period . All data come from the Central Institute of Statistics (ISTAT) database according to the report on equitable and sustainable well-being (BES). Given the hierarchical structural conformation of the model and a sufficiently large number of years describing socio-economic and health issues, it is able to capture: (i) endogeneity issues; (ii) interdependency, commonality, and homogeneity; (iii) heterogeneous effects; and (iv) model misspecification problems.

According to the TROB procedure, in the first stage, I find that 39 best covariates better fit the data with Posterior Inclusion Probabilities (PIPs) in (10) and . Thus, I obtain best model solutions (). Because of the curse of dimensionality, I further shrink the data performing the second stage involved in the TROB procedure. In Table 1, I display the final 22 best covariates, obtaining best model solutions (), with PIPs in (11). All variables and descriptive statistics are displayed in Tables 4 and 5 on Appendix A, respectively.

Table 1

Best Candidate Predictors – TROB (second stage)

Idx.	Predictor	Label	Unit	PIP(%)	CPS
Healthcare and social statistics
1	Overweight	Overobe	Std. rates per 100 people	44.13	0.002
2	Obesity	Obe	Std. rates per 100 people	45.74	0.002
3	Consumption of tobacco	Smoke	Std. rates per 100 people	37.15	0.007
4	Consumption of alcohol	Alcohol	Std. rates per 100 people	0.43	0.034
5	Sedentary rate	Sed	Std. rates per 100 people	11.51	0.019
6	Appropriate use of food	Food	Std. rates per 100 people	43.62	0.047
7	Cultural interests	Culti	% values	0.36	0.873
8	Social participation	Socpar	% values	21.05	0.996
9	Risk of poverty	Rop	% values	27.35	0.001
10	Public transport satisfaction	Ptran	% values	0.38	0.830
11	Health care expenditure	Hexp	Thousands € (log.)	31.53	0.782
Economic status
12	High school diploma	Hschool	% values	0.41	0.864
13	Graduates/other qualifications	Degree	% values	55.37	0.999
14	Employment rate	Employ	% values	43.24	0.975
15	Fixed-term contract (= 5 years)	Fterm	% values	49.15	0.002
16	Weighted income per capita	Income	Thousands € (log.)	0.38	0.947
17	expenditure on R &D	Exprd	Thousands € (log.)	24.33	0.986
Individual-specific factors
18	Neither studing nor working	Nsw	% values	31.28	0.008
19	Family relationship satisfaction	Family	% values	66.12	0.983
20	Free time satisfaction	Ftime	% values	37.50	0.991
21	Poorly-paid job	Ppjob	% values	36.74	0.001
22	Non-regular employment	Nnemploy	% values	0.17	0.022
–	Added value per employee	Ave	Thousands € (log.)	–	–

The Table is so split: the first column denotes the predictor number; the second and the third column describe the predictors and the corresponding labels; the fourth column refers to the measurement unit; and the last two columns display the PIPs (in ) and the CPS, respectively. The last row refers to the outcomes of interest. The contraction std. stands for ’standardized’. All data refer to ISTAT database

Table 4

Data source – subset

Idx.	Variable	Source
Healthcare and social statistics
1	Overweight	ISTAT: Aspects of Daily Life Survey
2	Obesity	ISTAT: Aspects of Daily Life Survey
3	Consumption of tobacco	ISTAT: Aspects of Daily Life Survey
4	Consumption of alcohol	ISTAT: Aspects of Daily Life Survey
5	Sedentary rate	ISTAT: Aspects of Daily Life Survey
6	Appropriate use of food	ISTAT: Aspects of Daily Life Survey
7	Cultural interests	ISTAT: Aspects of Daily Life Survey
8	Social participation	ISTAT: Aspects of Daily Life Survey
9	Risk of poverty	ISTAT: ’EU-SILC’ Survey
10	Public transport satisfaction	ISTAT: Aspects of Daily Life Survey
11	Health care expenditure	ISTAT: Aspects of Daily Life Survey
Economic status
12	High school diploma	ISTAT: Labour Force Survey
13	Graduates/other qualifications	ISTAT: Labour Force Survey
14	Employment rate	ISTAT: Labour Force Survey
15	Fixed-term contract (= 5 years)	ISTAT: Labour Force Survey
16	Weighted income per capita	ISTAT: National Accounts
17	Expenditure on R &D	R &D in Enterprises Survey
Individual-specific factors
18	Neither studing nor working	ISTAT: Labour Force Survey
19	Family relationship satisfaction	ISTAT: Aspects of Daily Life Survey
20	Free time satisfaction	ISTAT: Aspects of Daily Life Survey
21	Poorly-paid job	ISTAT: Labour Force Survey
22	Non-regular employment	ISTAT: National Accounting
-	Added value per employee	Small- and Medium-sized Enterprises (SMEs) Survey

Here, ’EU-SILC’ stands for ’European Union-Income and Living Conditions’ Survey (https://www.istat.it/en/archivio/216934’eu-silc.istat’). All data come from ISTAT database (https://www.istat.it/en/archive’documents_archive’)

Table 5

Descriptive statistics – subset

Idx.	Predictor	Mean	Standard deviation	Median	Skewness	Kurtosis
Healthcare and social statistics
1	Overobe	45.27	3.66	43.60	0.57	-1.49
2	Obe	55.07	3.64	53.50	0.52	-1.46
3	Smoke	21.46	1.63	21.60	0.20	-0.71
4	Alcohol	18.13	3.28	18.50	0.19	-0.90
5	Sed	41.18	10.29	39.60	0.26	-1.51
6	Food	19.10	3.58	20.90	-0.56	-1.47
7	Culti	28.76	6.45	31.60	-0.58	-1.40
8	Socpar	24.17	4.20	25.20	-0.37	-1.25
9	Rop	19.88	9.99	15.40	0.58	-1.57
10	Ptran	14.48	4.14	13.10	0.66	-0.80
11	Hexp	59.83	8.32	64.50	-0.71	-1.45
Economic status
12	Hschool	57.41	6.91	58.40	-0.28	-1.18
13	Degree	22.67	4.99	22.10	0.23	-1.03
14	Employ	61.33	9.82	66.00	-0.61	-1.49
15	Fterm	19.73	4.73	18.50	0.29	-1.42
16	Income	4.24	0.09	4.28	-0.56	-1.56
17	Exprd	1.20	0.25	1.30	-0.41	-1.43
Individual-specific factors
18	Nsw	21.97	8.01	19.60	0.46	-1.30
19	Family	33.87	5.13	33.60	0.05	-1.31
20	Ftime	64.40	4.06	65.80	-0.45	-1.19
21	Ppjob	11.45	4.28	9.90	0.41	-1.57
22	Nnemploy	13.60	3.88	12.60	0.34	-1.54
–	Ave	47.31	9.20	46.10	0.18	-1.16

Descriptive statistics for outcomes of interest and best subset of covariates involved in the analysis are shown. They account for the 21 Italian regions spanning the period . All data come from ISTAT database

Here, some considerations are in order. (i) Most model uncertainty and overfitting are avoided: indeed, the Conditional Posterior Sign (CPS) deals well with the sign certainty, taking values close to 1 (or 0) if a covariate in has a positive (or negative) effect on wage. (ii) Healthcare and social statistics matter as much as economic variables when evaluating cross-sectional wage effects (similar PIPs). (iii) Individual-specific factors and exogenous indicators need to be assessed when studying the effects of socio-economic factors on labour market outcomes (highly large PIPs). Cross-sectional individual factors are supposed to have homogeneous effects on wage. Contrarily, all non-homogeneous effects would be captured by via nonparametric Bayesian inference. (iv) The strictly exogenous factors () contain dummy variables to assess and then test the presence of structural breaks, evaluated through the Chow test, to investigate further interdependency among covariates and outcomes when performing density forecasts and dealing with (potential) spatial spillover effects. Indeed, labour market outcomes at a regional level are usually characterized by spatial dependence and then overweight may also exhibit spatial interdependence (see, for instance, [26] and [34]). (v) Finally, all predictors with PIPs displayed in bold correspond to the ones to be accounted for the final solution.

Nevertheless, although the intense shrinkage in the parameter space, the final solution would still require more effort: indeed, there are 22 potential best (endogenous) predictors better fitting the data. Thus, I run the last stage involved in TROB procedure to obtain the final model solution better performing the data with higher log Bayes Factor. It consists of 16 final best covariates with and so split: predictors (1, 2, 3, 5, 6, 8, 9, 11) for ; predictors (13, 14, 15, 17) for ; and predictors (18, 19, 20, 21) for . All their lags are then involved as instruments in the estimating procedure to deal with endogeneity issues and functional forms of misspecification, with particular emphasis on the lagged outcomes () and predictors 1, 2, 5, 11 addressing (potential) causal effects and spatial dependence, respectively (Table 2). In the dynamic panel data regression, I also include three time-invariant effects (, , and ), tested through the Chow test, denoting the presence of (potential) structural breaks in: 2007, due to the global financial crisis; 2011, due to fiscal consolidations in post-financial crisis; and 2014, due to recovery funds5. Moreover, because of possible significant difference between males and females affecting outcomes, I run three different models: (i) Model 1 accounting for all individuals i’s; (ii) Model 2 focusing on male individuals; and (iii) Model 3 accounting for female individuals. Gender disparities have been tested through F-test by comparing the variances of the two samples. More precisely, I test whether there are no differences among the two variances for males and females () versus the two-sided alternative (). Overall, the coefficients differ between males and females with: ; ; and ratio of explained variances equals to , where T stands for ’Total variance’.

Table 2

Dynamic regression results

Idx.	Variables	Model 1 - total	Model 2 - male	Model 3 - female
Healthcare and social statistics
1	Overobe	-1.84∗∗(0.83)	-0.98∗∗(0.43)	-1.23∗∗(0.54)
2	Obe	-1.92∗∗∗(0.91)	-1.16∗∗∗(0.51)	-1.21∗∗∗(0.52)
3	Smoke	-1.33∗∗(0.52)	-0.94∗(0.46)	-1.35∗∗(0.53)
5	Sed	-0.57∗∗(0.23)	-0.68∗∗(0.25)	-0.81∗∗∗(0.21)
6	Food	-1.25∗∗∗(0.16)	-1.15∗∗∗(0.21)	-1.32∗∗∗(0.18)
8	Socpar	0.52∗(0.30)	0.32∗∗(0.16)	0.25∗∗(0.15)
9	Rop	-0.57∗∗(0.26)	-0.44∗(0.25)	-0.46∗∗(0.28)
11	Hexp	1.35∗∗∗(0.65)	1.46∗∗∗(0.73)	1.49∗∗∗(0.75)
Economic status
13	Degree	0.92∗∗∗(0.30)	1.11∗∗∗(0.29)	0.61∗∗∗(0.21)
14	Employ	3.77∗∗∗(0.26)	3.77∗∗∗(0.33)	3.72∗∗∗(0.29)
15	Fterm	-1.94∗∗∗(0.32)	-1.68∗∗∗(0.25)	-1.74∗∗∗(0.40)
17	Exprd	9.56∗∗∗(2.61)	8.91∗∗∗(3.22)	8.53∗∗∗(2.99)
Individual-specific factors
18	Nsw	-2.42∗∗∗(0.34)	-2.26∗∗∗(0.44)	-2.31∗∗∗(0.38)
19	Family	0.78∗∗∗(0.24)	0.73∗∗(0.31)	0.75∗∗(0.35)
20	Ftime	0.44∗∗∗(0.13)	0.37∗∗(0.16)	0.37∗∗∗(0.13)
21	Ppjob	-2.75∗∗∗(0.53)	-2.69∗∗∗(0.47)	-2.72∗∗∗(0.52)
Exogenous factors
–	x_1t	-1.09∗∗∗(0.25)	-1.14∗∗∗(0.026)	-1.33∗∗∗(0.30)
–	x_2t	1.55∗∗∗(0.38)	1.66∗∗∗(0.38)	1.58∗∗∗(0.37)
–	x_3t	2.18∗∗∗(0.42)	2.26∗∗∗(0.44)	2.17∗∗∗(0.43)
Lagged predictors - spatial dependence
–	L. overobe	-1.89∗∗(0.62)	-1.74∗∗∗(0.54)	-1.67∗∗∗(0.55)
–	L.obe	-2.03∗∗∗(1.12)	-1.75∗∗∗(1.09)	-1.79∗∗∗(1.11)
–	L.sed	-0.67∗∗(0.34)	-0.73∗∗(0.38)	-0.76∗∗∗(0.40)
–	L.hexp	1.59∗∗∗(0.82)	1.63∗∗∗(0.83)	1.65∗∗∗(0.84)
Lagged outcomes - dependent variable
–	L.ave	0.54∗∗∗(0.07)	0.51∗∗∗(0.09)	0.55∗∗∗(0.08)

The Table is so split: the first two columns denote the predictor number and labels; and the last three columns display the estimates and the standard errors (in parenthesis). Here, the variable of interest is the added value per employee and L. stands for the first order lag operator concerning the outcomes and some endogenous predictors to deal with causal relationships and spatial dependence, respectively. The significant codes are: *** significance at 1; ** significance at 5; and * significance at 10

The DPD estimates are displayed in Table 2. Concerning Model 1, six findings are addressed. (i) Health status and social factors strongly affect wage rates and are essential in determining a stable and well-functioning economy, high-quality jobs, and then the stability of a country. The negative relationship between health sector and wage (predictors 1, 2, 3, 5, 9) highlights the positive impact that the former has on the economic performance in the national economy through the jobs it generates and the purchase of goods and services. (ii) Healthcare factors would reduce social exclusion at the local level (significance and positive magnitude for predictor 8) due to their impact on employment, working conditions, and household income. Thus, the health sector can increase economic status acting as a key sector for driving forward the implementation of local and national goals for sustainable development. For instance, health and social work activities have constituted around of total employment in 2015 according to the Organisation for Economic Co-operation and Development (OECD) countries. (iii) Accounting for economic and individual-specific factors, the availability of good jobs lies at the heart of inclusive and sustainable growth eliminating poverty. More precisely, good working conditions would provide: decent pay ensuring at least a minimum or living wage (positive magnitude for predictors 19, 20); employee benefits, such as maternity and paternity leave, ensuring part- and full-time workers to receive similar benefits (negative magnitude for predictor 21); minimal use of temporary contracts (negative magnitude for predictor 15); safe working conditions (negative magnitude for predictor 18); and opportunities for progression and career development (positive magnitude for predictors 13, 14, 17). (iv) Overweight (predictor 1), obesity (predictor 2), and sedentary rate (predictor 5) are likely simultaneously determined in the estimating procedure to deal with spatial dependence and produce unbiased, consistent, and efficient estimators. Here, simultaneity arises when a predictor is endogenous to the system and then is likely correlated with the error term. (v) According to the previous point, let the PIPs and the degree of relationship between wage effects and excess weight (predictors 1 and 2) be very close, it highlights how environmental and policy issues, including interventions to change the local environment to create opportunities for physical activity, would positively affect the economic dynamics and then its performance in the national economy. Moreover, the of obese individuals would tend to matter slightly more than the of overweight ones, emphasizing the negative impact that an excessive weight has not only on work activities but also on the health system (e.g., higher magnitude of excess lifetime medical care costs associated with an increasing degree of obesity). These findings find confirmation with related previous studies such as [70] among others. (vi) Finally, concerning time-fixed effects, policy programs () positively affect outcomes with larger magnitude than post-crisis fiscal measures (). Conversely, negative effects on outcomes are displayed during the Great Recession ().

According to gender (Models 2 and 3), three main results are in oder. (i) First, employment penalties persist for males and even more for females. The undervaluing of women’s social care and health jobs relate to wider norms and attitudes in society, where women’s skills are often less visible than men’s ones and harder to be quantified, mainly when linked to productivity. Generally, women are stereotyped as being naturally good at the job and are thought to be prepared to trade lower pay for job reward. Despite the continuing reduction in the disparity between male and female employment rates, large gender differences still remain in Italy. These findings recognize the results found in previous studies such as [11, 24], and [30]. (ii) Second, healthcare and social factors negatively affect employment opportunities and wage rates for males and slightly more for females. For instance, negative impacts on working conditions, e.g., due to smoking, alcohol use, and unappropriate use of food, would matter more for females. (iii) Third, another important factor to be accounted for understanding the pay levels of health and social care workers is that many more women than men work part-time hours in employment. Indeed, individual-specific differences highlight that discrimination still accounts for a significant amount of the gender pay gap, referring to the relative differences in average gross hourly earnings of women and men. However, the gender pay gap needs to be evaluated in either the economic status or healthcare and social context to understand the position of women in society. More precisely, in specific occupation rates of pay, women are likely to work in jobs that undervalue their skills and training, and to be numerically lower than in those where men predominate. These results are consistent according to the structural indicators drawn up by the European Union to monitor the implementation of the European Strategy for Jobs and Growth (see, for instance, the ’European Commission in the Roadmap for Gender Equality, ’). According to time-fixed effects, the results are similar to the ones found in Model 1, with a slightly larger (lower) magnitude for females according to negative (positive) effects on outcomes.

In Table 3, I display the main diagnostic tests dealing with significance and robustness of the estimating procedure, autocorrelation problems, panel unit root tests, spatial dependence, and structural breaks. Here, some considerations are in order. (i) The dynamic panel regression shows a and then is robust evaluating most of the explained variability of the outcomes of interest, with referring to the adjusted R2. (ii) the Sargan-Hansen’s test for over-identification () highlights the performance and usefulness of the DPD to deal with endogeneity issues and functional forms of misspecification ( [68] and [37]). Moreover, it is also able to address variable selection problems avoiding model uncertainty and overfitting (for instance, the CPSs close to 1 or 0 in Table 1). (iii) There are no serial correlations in idiosyncratic errors according to the first order lag operator. The test statistic refers to [6] (). (iv) The Multivariate Ljung-Box test () shows no serial correlation among residuals over time and then the time-series results are robust and valid ( [49]). (v) The highly large and Sargan’s test statistics would also highlight the importance to account for heterogeneous individual characteristics and socio-economic factors when studying the effects of healthcare status, overweightness, and labour market outcomes. These relationships are emphasized by testing the spatial dependence through the [5]’s robust Lagrange Multiplier (LM) statistics on the autoregressive coefficients and in (1). The statistics are distributed as a with one degree of freedom, where the null stands for no spatial autocorrelation. The results are robust to the presence of spatially lagged correlation in the variable of interest (ave) when testing for spatial correlation in excess weight (overobe and obe), and healthcare statistics (sed and hexp). (vi) Concerning panel unit root tests, I use three different test statistics for identifying unit root processes and trend-stationarity for every AR(1) time series involved in the system as instrument: (1) [35] () that extends the [43]’s test to panel data with i.i.d. heteroskedastic and serially correlated disturbance terms across i and over t by proposing a residual-based Lagrange multiplier test; (2) [39] () that suggests standardized t-bar test statistic based on the Augmented Dickey-Fuller statistics averaged across units; and (3) [52] () that proposes a comparative study of unit root tests with panel data. In the first study, stationarity around a deterministic trend (or level) is the null hypothesis against the alternative of a unit root. Conversely, in the last two studies, the unit root is the null hypothesis to be tested against the alternative of trend-stationarity. In all these test statistics, either the critical values or the moments of the asymptotic distributions have been evaluated using MCMC algorithms. The results confirm no unit root in the lagged time series and then likely ensure trend-stationarity. (vii) Finally, (potential) structural breaks matter for every estimation model affecting the outcomes. They are dealt with the three time-fixed effects included in the estimating procedure through F-statistics ().

Table 3

Diagnostic panel tests

Test statistic	Model 1 - total	Model 2 - male	Model 3 - female
Significance and Robustness
R¯^2	0.89	0.85	0.86
Q_SH	255.03 (0.00)	157.02 (0.00)	197.01 (0.00)
Q_A	1.78 (0.19)	2.86 (0.11)	0.99 (0.32)
Q_LB	0.02 (0.89)	0.07 (0.79)	0.04 (0.84)
Spatial dependence
Robust LM on overobe	9.97 (0.00)	11.30 (0.00)	11.64 (0.00)
Robust LM on obe	10.14 (0.00)	12.07 (0.00)	12.74 (0.00)
Robust LM on sed	5.44 (0.03)	5.83 (0.02)	6.03 (0.00)
Robust LM on hexp	12.87 (0.00)	13.12 (0.00)	13.46 (0.00)
Panel unit root tests
U_H	-4.16(0.13)	-4.35(0.21)	-4.68(0.23)
U_I	-3.10(0.09)	-4.16(0.02)	-3.71(0.03)
U_M	-3.38(0.06)	-3.23(0.08)	-3.42(0.07)
Time-fixed effects
F_x1	31.43 (0.00)	22.86 (0.00)	25.53 (0.00)
F_x2	16.04 (0.00)	17.73 (0.00)	15.59 (0.00)
F_x3	15.22 (0.00)	4.92 (0.00)	5.74 (0.00)

The Table refers to the results of diagnostic tests for significance and robustness in the estimating procedure, autocorrelation problems, panel unit root tests, spatial dependence, and structural breaks

Density forecasts and policy-relevant strategies

In Fig. 1, I draw conditional density forecasts through accurate MCMC algorithms and implementations. More precisely, posterior distributions are computed according to Eqs. (28)–(30) for and given , (34)–(35) for moment distributions in given initial values (, ), (31)–(32) for the final best parameter and model space, and (33) for . All data are expressed in standard deviations, and the accuracy is measured by Eq. (36).

Fig. 1

The plot draws density forecast combinations for wages across 21 Italian regions, spanning the period 2006 to 2021. It accounts for the final submodel solution . The Y and X axes represent the conditional projections and sampling distribution in years, respectively

Density forecasts are performed running 100, 000 iterations per each random start, spanning the period 2006 to 2021. The optimal number of draws has been chosen according to the lBF: indeed, with 100 and 1000 iterations, I obtain highly lower lBF (4.89 and 5.34, respectively). The associated computational costs are minimized, ensuring consistent posterior estimates and dimension reduction. Density forecast combinations correspond to the projections of every subsequence drawn in the sample. In this context, I divide the sample into three groups to better highlight the strengths and limits of policy strategies: (1) Northern Italy; (2) Central Italy; and (3) Southern Italy. The yellow and red curves denote the confidence bands, and the blue and purple curves denote the conditional6 and unconditional7 projection of outcomes for each time period t, respectively. Here, the outcomes absorb the conditional forecasts computed for a time frame of 3 years in order to also address potential findings concerning the impact of the ongoing pandemic crisis on the global economy. The natural conjugate prior refers to the three subsamples according to the structural breaks involved in the panel regression: 2007; 2011; and 2014.

The results emphasize the findings obtained in Table 2. (i) Conditional projections lie in the confidence interval; conversely, unconditional projections tend to diverge over time. Thus, when studying time-varying healthcare and socio-economic factors, structural time-fixed effects () and correlated individual-random effects () have to be accounted for. (ii) Density forecasts in group 1 (top-plot) tend to show a similar evolution, but with a larger spillover effects to the ones in group 2 (middle-plot) because of not directly observed regional and European country interdependencies (e.g., similar dynamics, co-movements, economic structure). (iii) Density forecasts in group 3 (bottom-plot) tend to display larger divergences, in terms of spillovers, because of potential functional forms of misspecification (e.g., strong inter-linkages across regions in Northern Italy driving the transmission of unexpected shocks on labour market outcomes). (iv) Most positive projections in groups 1 and 2 higlight the accuracy of the methodology by dealing with endogeneity issues. For instance, regions in the Northern Italy tend to be net sender because of consistent interdependencies among them and Northern Europe, and then Central Italy regions because of stringent economic-institutional linkages with the former. (v) The spillover effects show a sudden change in their dynamics according to the three structural breaks evaluated in the panel data regression. Thus, it highlights the presence of persistent (unobserved) shocks over time affecting time-varying predictors and then outcomes. (vi) Accounting for the ongoing pandemic crisis (conditional projections in ), a sufficiently stronger recovery seems to matter across the regions in Northern Italy. Conversely, slower socio-economic improvements would affect the regions in Central Italy and, even more, the southwards. Thus, the need for forecasters and policymakers to deal with cross-sectional individual-specific heterogeneity when investigating dynamic feedback between socio-economic effects and labour market outcomes.

From a policy perspective, the above empirical results can be so summarized. (i) Appropriate prevention policy strategies in terms of higher healthcare and social statistics need to be improved, mainly in accordance with the current pandemic disease. Indeed, solid and coordinated improvements in healthier and safer workplace would increase employment opportunities and wage effects, achieving higher labour productivity and (possibly) lessening social disparities in gender. (ii) Better and more active workforce should be designed to deal with key behavioural and individual-specific risk factors, such as obesity, social participation, non-regular employments, and leisure satisfaction, negatively affecting employment opportunities and wage rates for males and slightly more for females. (iii) Since 2011, a substantial number of labour market and social policy reforms has been achieved, improving the flexibility of job creation and extending social assistance and social protection coverage. In 2014, further measures have been enacted under the so-called Jobs Act reform (’Renzi government’). The main goal of the Jobs Act reform has been twofold. First, it has aimed to encourage jobs growth by making labour input cheaper. Second, the reform has aimed to make contracts more appealing to employers by the introduction of ’increasing-protection employment contracts’. These refer to new permanent open-ended contracts, with lower dismissal costs and reduced legal uncertainty compared to the previous ones to promote a combination of labour market flexibility and security for workers. (iv) Recently, in 2018, a partial reversal of the Jobs Act provisions has been carried out by the Conte government, the so-called ’Dignity Decree’, by restricting the use of fixed-term contracts and increasing the costs of dismissals. According to social policies, a combination of vocational education and employment measures has been also revised under the so-called ’Good School reform’, by including initiatives to simplify short-term contracts and then align them better with labour market needs. Nevertheless, besides the weak labour market recovery in terms of job creation, the Italian labour market is still characterised by increasing instability of employment and use of fixed-term contracts.

From a global perspective, related considerations are in order. (i) European Union Structural Funds have supported the setting up the necessary infrastructure in terms of services and operators. More precisely, the ’Stability Law’ in 2016 have introduced the ’Support for Active Inclusion’ (SIA) to address the way for the new ’Inclusion Income’ implemented in December 2017. However, due to the lack of resources, the new measure has been not sufficient to cover all individuals living in absolute poverty, and the monetary transfer to families follows to be still very small. (ii) Concerning the COVID-19 pandemic crisis, confirmed in Italy on 31 January 2020, it has rapidly expanded across regions in Northern Italy and then the southwards, although the impact at the local level has been heterogeneous (see, for instance, Figure 1). On 22 April 2020, Italy has been the third country in the world by the number of reported COVID-19 infection cases, and the second by the number of deaths among the infected patients (see, for instance, the ’Center for Systems Science and Engineering’). The progressive tightening of containment measures has affected more and more larger segments of the Italian economy, mainly with regards to retail trade businesses with the exception of those related to the sale of food and other basic necessities. Since the beginning of March 2020, the OECD has designed relevant policy measures aimed at improving the healthcare system and providing economic support to households, workers, and businesses. Then, the second set of strategies has been also developed supporting employees, self-employed, and businesses. (iii) Finally, according to the European Union’s regional policies (2018) and the post-2020 Cohesion policy reform, several ’place-based’ programs have been addressed for the redistribution of wealth across regions and countries to stimulate growth in the areas lagging in development. In that context, three main objectives can be summarized. First, public finances and current expenditure should be controlled but allowing development policies to be implemented in backward areas with the greatest potential. Second, an adequate system of fiscal compensation has to be included to offset the structural competitive disadvantages to which the Southern Italy regions are exposed, mainly within the Eurozone, with a view to progressively overcoming them. Third, a rebalancing of the current geopolitical configuration to aim for cooperation and development policies for the Southern area needs to be accounted for.

Concluding remarks

This study improves the existing literature when investigating the effects of obesity, socio-economic variables, and individual-specific factors on work productivity across Italian regions. A dynamic panel data with correlated random effects is used to address variable selection problems in high dimensional time-varying data. Semiparametric Bayesian inference is involved in the estimating procedure to jointly deal with incidental parameters, endogeneity issues, and functional forms of misspecification. The main thrust of the proposed methodology is the use of conjugate hierarchical informative (proper) priors to discover the most probable subset of potential predictors affecting outcomes to capture different and/or not directly observed dynamics and interconnections among lagged endogenous factors. Posterior distributions for conditional density forecasts are obtained by MCMC algorithms and implementations. Cross-sectional Heterogeneity is modelled nonparametrically allowing for correlation between heterogeneous parameters and initial conditions as well as individual-specific regressors.

An empirical example describes the estimating procedure and forecasting performance. More precisely, the estimation sample accounts for all the 21 Italian regions covering the period . All data come from the Central Institute of Statistics database according to the report on equitable and sustainable well-being. Conditional density forecasts are then performed spanning the period 2006 to 2021 to also address potential findings concerning the impact of the ongoing pandemic crisis on the global economy. In that context, the sample is divided in three groups to better highlight the strengths and limits of policy strategies: (1) Northern Italy; (2) Central Italy; and (3) Southern Italy. Moreover, natural conjugate priors are involved in the procedure to also deal with three permanent time-invariant effects affecting panel data: 2007, due to the global financial crisis; 2011, due to fiscal consolidations in post-financial crisis; and 2014, due to recovery fund.