Modeling unemployment duration, determinants and insurance premium pricing of Malaysia: insights from an upper middle-income developing country

This paper investigates the behavior and determinants of unemployment in Malaysia, an upper middle-income country with relatively low unemployment but rising labor market mismatches, and derives a model-based insurance premium pricing schedule that broadens the country's unemployment insurance coverage. The unemployment profile is best fitted with Burr distribution, and using survival analysis, the average unemployment duration is found to be the highest for the youths and it decreases gradually after age 20. The results of the Cox proportional hazards model show that age, gender, and education are significant explanatory variables of unemployment. While all age groups are statistically significant in explaining unemployment, the youngest cohort has the least likelihood of being hired, while the older groups are found to have a lower hazard ratio that could be attributed to work experience. The significant gender effect is also characterized by the female workforce taking a longer time to exit unemployment as shown by its lower hazard rate than male workers. Similar to findings for other countries, those with higher education have the highest chance of gaining employment. Fitting the estimated parameters to a standard unemployment insurance premium pricing framework based on the Capital Asset Pricing Model (CAPM), the analysis shows that the coverage of the existing national unemployment insurance scheme can be expanded to all income categories at lower premium prices. The widening of the coverage will contribute to the strengthening of the country’s social safety net in addition to reducing inequalities in the post-COVID-19 period. Supplementary Information The online version contains supplementary material available at 10.1007/s43546-022-00287-2.

Introduction
Unemployment remains a primary and recurring concern in most countries regardless of development stage and unemployment level. The multifaceted issues range from massive job losses and anemic recovery following the 2007-2009 global financial crisis to the high youth unemployment in the middle-eastern countries that has been cited as a key factor for the Arab Spring uprising. More contemporary issues include unequal access to job opportunities in low unemployment economies and the potential large job displacements caused by automation and digital technologies and skills mismatches in middle and high-income countries.
An upper-middle-income country since the early 1990s, Malaysia's unemployment rate has been relatively low and stable, averaging 3.3% over the last 2 decades. By contrast, the average rate for upper-middle countries is 5.6%, while that for highincome nations is higher at 6.8%. Notwithstanding the country's low unemployment rate, the Malaysian government introduced national unemployment insurance in 2017. There is however limited statistical modeling of unemployment behavior both in Malaysia and the developing world.
This study fills the gap by estimating Malaysia's unemployment distribution and duration, two crucial inputs in the design of unemployment insurance schemes. The study also examines key factors influencing unemployment. Despite the country's low unemployment rate, a comparison of key determinants with those of countries with high unemployment rates will have significant policy implications in addressing unemployment issues that are common as well as unique to each country.
The objectives of the unemployment analysis are twofold. First, this paper provides a deeper insight into the patterns of the unemployment rate. The data analysis shows that the Burr distribution form is appropriate for the model fitting as it outperformed other positively skewed distribution functions. Second, the survival analysis via hazard function corresponding to Burr distribution is used to estimate unemployment duration which is a key input in pricing unemployment insurance premium. To motivate the enhancement of the country's recently launched national unemployment insurance scheme, the study applies the derived statistical estimates to estimate a model-based unemployment insurance premium rates by wage classes.
In survival analysis, the hazard function defines instantaneous death. We apply this definition to explain unemployment duration. Unemployment duration is defined by the conditional probability of exiting unemployment, which means that an unemployed person will find employment at time t, given that at time t-1, the person was unemployed. This is illustrated by the hazard function or commonly known as hazard rate. This paper explores the relationship between two covariates. The hazard ratio is used to measure the relativity. If the ratio is equal to 1, there is an equal possibility of being unemployed at time t-1 and employed at time t. If it is more than 1, it means that the employment at time t has a higher chance to happen, while for a ratio lower than 1, the situation of being unemployed will remain longer than being employed. This statistical approach is employed to analyze the underlying unemployment situation in Malaysia and provide an empirical estimation of unemployment duration, determinants, and insurance premium pricing to support evidence-based policy decisions.
The rest of the paper is organized as follows. The relevant literature is reviewed in the section "Literature review". The section "Data sources and description" describes the labor survey data sourced from the Department of Statistics Malaysia (DOSM). The section "Modeling and estimation methods" details the model and estimation methods for fitting unemployment rates and deriving unemployment duration using the Burr distribution function. This section also explores covariates relationships using Cox proportional hazard regression models and enumerates the construction of the unemployment insurance scheme based on the Equivalence Principle. The section "Results and discussion" discusses the estimation results and the section "Conclusion" concludes.

Literature review
Unemployment is intertwined with a country's macroeconomic management where it is commonly aimed at securing sustained growth, price stability, and full employment. Unemployment below the natural rate reflects underutilization of a country's human resources and, consequently, lower gross domestic product (GDP) growth potential. Besides a loss of tax revenue to the government, unemployment is also associated with social and political problems such as crime, vandalism, and popular uprisings. Among the earliest to relate unemployment to macroeconomic impact is Malinvaud (1982) who suggested a national unemployment insurance system to minimize unemployment risk. Beenstock (1985) enhanced the unemployment insurance scheme to further consider an individual's characteristics, specifically, working behavior. He proposed that the premium paid should fairly reflect different risk groups. He then constructed a competitive pricing model for unemployment insurance. The model was subsequently enhanced with stochastic properties by Blake and Beenstock (1988), taking inflation into the consideration as a factor in determining premium pricing. Other studies such as Pollak (2012) established the relationship between human capital depreciation and unemployment insurance policy.
Some existing study is found applying survival analysis to analyze the unemployment spells and the measurement of the effects of the unemployment condition. Lakuma et al. (2016) applied survival analysis for a regional unemployment study in Uganda. They reported that the cohort with higher education surprisingly gains employment slower than those with lower education across the region. Women have a lower possibility to exit unemployment and the chances of being employed increase according to age. While Malaysia is advisable to focus on skill development for the age cohort of 19-21, in Uganda, it is suggested that the investment in training skills should be given to the cohort younger than 25 years old. Bajram (2013) found that re-employment happens in the first 59 months with the frailty model in survival analysis. Bowers (1980) applied transition probabilities to estimate the unemployment duration. In the transition probabilistic model, he defined three states-employed, unemployed, and not in the labor force. Salant (1977) discussed the Pareto distribution, while Chuang and Yu (2010) proposed Weibull distribution for estimating the unemployment duration in Taiwan's labor force. In an enhancement to the estimation technique, Simwa et al. (2016) applied a mixture of the Weibull distribution to accommodate the heterogeneity in the US unemployment data. The non-parametric Kaplan-Meier estimator has been used to estimate the unemployment duration.
Given that the existing literature is limited and focused primarily on developed economies, this paper seeks to address the lack of statistical modeling in unemployment studies for developing countries. This research provides an opportunity to extend the coverage using Malaysia as a case study. The country's status as an uppermiddle-income country with a labor market characterized by low unemployment but a high number of unskilled foreign workers could shed light on the structural characteristics influencing unemployment behavior. As unemployment is an uncertain event, one can consider applying probabilistic models to measure it numerically. Unemployment duration, defined as the length of time in which a person seeking a job remains unemployed, therefore can be modeled by conditional probabilities. We applied a parametric approach in this study where the Burr distribution is fitted to the data set given its superiority over other forms in terms of flexibility.
While there are limited probabilistic model-based studies on unemployment, several researchers have applied other estimation approaches. Aitkin and Healey (1985) applied logit models to unemployment rates. Pereira et al. (2018) focused on the application of Bayesian hierarchical models to the Portuguese labor force survey data and made comparisons with the multinomial model and Beta model. Time-series analysis and forecasting are popular approaches to modeling unemployment. Holmes et al. (2013) applied a high-dimensional vector autoregressive analysis to study the behavior of unemployment rates in the US over time and across space. Christos (2005) forecasted the UK unemployment rate with 30 years data set. In time-series modeling, the researchers found that the autoregressive conditional heteroscedasticity models yield the best forecasts. Claveria (2019) forecasts the unemployment rate in eight European countries with autoregressive integrated moving average models. Kreiner and Duca (2020) show that the machine learning approaches outperform time-series models to provide better forecasts at shorter horizons. When applied to actual data, his results indicate that the machine learning model can identify the turning points of the unemployment rate. In comparison, the survival analysis remains the advantage to analyze the unemployment issue as it can measure the effect of the covariates.
Given Malaysia's low unemployment rate, a comparison of unemployment duration and the key determinants of unemployment with those of countries with high unemployment rates will have significant policy implications in addressing both unemployment issues that are common as well as unique to each country. Most researchers focus on the analysis of the determinants such as age, gender, education attainment level, and region. The east European countries, such as Romania, Bulgaria, Hungary, and Slovenia, are of interest, because they have a high unemployment rate among the countries in the region. Romania has an average unemployment rate of 6.64% from 1999 to 2018, among the lowest in the region. Bulgaria on the other hand recorded the highest average unemployment rate at 10.94% for the past 20 years. Hungary and Slovenia have recorded 7.36% and 6.95%, respectively. Danacica and Paliu-Popa (2017) investigated the determinants of unemployment in the East European countries where unemployment reached its highest level in 2011 following a recession. Focusing on explanatory variables such as gender and education, the authors found that in Romania, women experienced a 14% lower exit rate than men. In other words, it is more challenging for women to leave their unemployment state. The study also found that 22.8% of the unemployed are educated at the high school level, 9.1% with higher education, and 10.8% with unknown education. Those with a college education have the highest chance to find employment, followed in declining order by those with university education, high school, and unknown education. Similar results have been found by Ciuca and Matei (2011) where those with higher education have the lowest probability of unemployment. The study also revealed the 36-45 age cohort faced the most difficulty in finding employment, while it is easier for those aged 15-25 years. The study revealed that among the key variables tested, education has the highest power in explaining unemployment. In contrast, our finding for Malaysia shows education and gender as having an equal impact on unemployment with both groups facing the same possibility to exit unemployment. Kavler et al. (2009) reported the unemployment exit rate is 32% higher for men than for women in Croatia followed by Austria (29.9%), Slovenia (20.8%), and Romania (16.3%). Malaysia shares similar features whereby women take longer to exit unemployment and face more hurdles in finding jobs. The study on East European countries found older workers to be disadvantaged in finding employment in contrast to Malaysia where the unemployment duration is less than 30 days for the age group greater than 50 years. Another difference found in our study is that in Malaysia the unemployed with higher academic qualifications find it easier to leave unemployment as evidenced by the shorter unemployment duration for those with master's degrees compared to those with bachelor's degrees.
There are various other studies on the influence of factors such as age, gender, race, living region, and education attainment. In some other European countries. D' Agostino (2000) found that gender has an impact on unemployment in Belgium, Greece, France, Spain, Denmark, and Portugal, but not Italy. In Italy, gender does not exert much influence on unemployment. The effect of location or place of living area is prominent in Italy where those living in the south have greater difficulties in finding jobs. The study also shows the age variable has the least impact in Italy, the United Kingdom (UK), and Spain, but in Portugal, France, and Denmark the older cohorts find more difficulties in exiting unemployment. Higher education levels shorten the unemployment duration in the UK, Belgium, and Ireland, while in Greece and Spain, the educational level does not have a strong effect on the expected duration. On the other hand, Baussola and Mussida (2017) and Mussida and Fabrizi (2014) reported contradictory findings. They found gender has an impact on unemployment outflow.
In the UK, the analysis of the British Household Panel Survey by Boheim and Taylor (2000) indicates that female exit rates are longer than men. Women need 9 months to leave unemployment, whereas the men take only 6 months. Another key finding is that men aged over 25 are entering unemployment with an increased education level. Highly educated women are significantly more likely to enter part-time employment from unemployment than men. Most studies report that women face greater challenges in looking for a job, the exception being the study by Hoffman (1991) who found white females in the US are more likely to leave unemployment, particularly the higher educated white females who have better opportunities to secure jobs. Mirroring the finding for European countries, the US study found that older persons are less likely to exit unemployment.
There is a paucity of similar studies for developing countries where unemployment is expected to vary widely in scope and intensity depending on an individual country's economic performance, population growth, and labor market characteristics. Slow-growing economies face higher unemployment, while the more dynamic ones are likely to face different unemployment challenges.

Data sources and description
The data are sourced from the Labor Force Surveys conducted by the Department of Statistics once every 3 years for the period 1997-2015. There are 15 variables enumerated according to the International Standard Classification of Education (ISCED) classification that is analyzed, covering a total of 427,714 sample points. Table 1 provides the labor force count and participation rate by gender. The key variables are highlighted below: Age: A key variable in the employability modeling and analysis, the working age is based on the universally adopted range of 15-64 years, while youths are defined as those in the 15-24 age group.
Employment status: There are three categories in the labor force summary (see Supplementary Table 10 for the definition). The three labor force categories are Employed (code 1), Unemployed (code 2), and Outside Labor Force (code 3).
State: There are 13 states and 3 separately coded federal territories resulting in a total of 16 localities (see Supplementary Table 11 for the state label).
Educational attainment: There are four categories in the educational attainment variable as defined in Supplementary Table 12. The education system is divided into preschool education, primary, secondary, pre-university, and tertiary education. To enable a more focused analysis this study employs four categories: no formal education (code 1), primary (code 2), secondary (code 3), and tertiary (code 4) education.
Gender: An equal proportion of males (code 1) and females (code 2) are drawn from the survey. Table 1 shows the labor force participation by gender. Highest certificate obtained: There are nine categories under this explanatory variable, comprising, primary school achievement test (UPSR) or equivalent (code 1), lower secondary evaluation (PT3/PMR/SRP) or equivalent (code 2), the Malaysian certificate examination (SPM) or equivalent (code 3), Malaysian higher school certificate (STPM) or equivalent (code 4), certificate less than 6 months (code 5), diploma (code 6), degree (code 7), no certificate (code 8), and not applicable (code 9).
Field of study: This variable consists of nine fields of study. They are general programs (code 1), education (code 2), humanities and arts (code 3), social sciences, business and law (code 4), science, mathematics and computing (code 5), engineering, manufacturing and construction (code 6), agriculture and veterinary (code 7), health and welfare (code 8) and services (code 9). Figure 1 shows the unemployment rate for the cohort with tertiary education.
The data have been filtered and cleaned before the analysis. The descriptive statistics on the mean, median, first quartile, third quartile, range, standard deviation, and skewness of the unemployment rate for each variable are shown in Table 2. The standard deviations are small. It means that the data are not likely to disperse away from the mean and so it is reliable. The findings show that the explanatory variable has an impact on unemployment. The effects of these explanatory variables are tested with the Cox proportional hazards model.

Modeling and estimation methods
This section describes the selection of the distribution specification for the unemployment rate followed by the application of hazard function and hazard ratio in survival analysis to estimate unemployment duration. To investigate the key factors influencing unemployment, the Cox proportional hazards modeling technique is used. The Equivalence Principle method is used to determine the premium in the insurance field based on the principle that the expected value of the future loss function is equal to zero.

Distribution specification for unemployment rate
It is observed that the distribution of Malaysia's unemployment rate is positively skewed, as shown in Fig. 2. Given that the Burr distribution is characterized by unimodality and strong positive skewness, it is tested against other unemployment distribution functions for the working-age population (15 -65 years). The probability density function (pdf) is given by where k > 0, > 0, c > 0 , k, , c are the parameters of Burr distribution. The parameters are derived using maximum-likelihood estimation. The random variable x is age considered throughout the year (time) from 1997 to 2015.

Hazard function for unemployment duration
Unemployment duration is estimated using the hazard function given by the mathematical expression In general terms, Eq. (1) gives the failure rate of an item given that the item has survived for a period of time t . In this study, the term 'survive' refers to being unemployed. Therefore, Eq. (1) is interpreted as the chance of an unemployed individual at x the age to be unemployed for another interval dt given that the person has been unemployed for a period of time t . This is the underlying idea of unemployment duration estimation.
To estimate the unemployment duration for a person aged x , we consider the alternative expression of the hazard function. For simplification, Eq. (1) is expressed as where f (t) is the probability density function for unemployment rate and S(t) is the survival function which represents no change in status, in other words, remaining unemployed.

Cox proportional hazards model
The Cox proportional hazards model was introduced by a British statistician, David Cox in 1972 when he published the paper entitled "Regression analysis and life tables". Since then, the use of the model has been extremely significant in biomedical research. Cox's proportional hazards model is mainly concerned with investigating the effects of the explanatory variables. The hazard function to measure the unemployment scenario is represented by where x i are covariates ( x 1 = age , x 2 = gender , x 3 = education , x 4 = ethnic , x 5 = state , x 6 = highest certificate obtained, and x 7 = field of study ), i (t) is the hazard function of the individual explanatory variables i , and 0 (t) is the baseline function which corresponds to an observation (unemployment rate) with no covariates. The quantiles exp i represents the hazard ratio. The hypotheses to test the significance of the effects of the respective coefficients are as follows: A critical p value of less than = 0.05 is adopted to reject the null hypothesis which implies that the coefficient is an influencing factor for unemployment.
The Hazard Rate (HR) is also used to measure the relationship between two covariates. An HR larger than one indicates that a covariate is positively associated with unemployment, an HR less than one indicates a covariate is negatively associated with unemployment, and an HR equal to one indicates no effect between the covariates and unemployment. In short, the HR can be interpreted by

Unemployment insurance premium pricing framework
The Equivalence Principle is a method to determine the premium in the insurance field. The idea behind the principle is that the expected value of the future loss function is equal to zero. It is applied here to determine the premium rate for unemployment insurance. An insurance policy is a financial agreement between an insurance company and its policyholders where the insurer provides security to the policyholder by paying some benefits during the event of an occurrence, as well as covering expenses associated with maintaining the policy contract. On the other hand, the insured consents to pay a premium to the insurance company to secure these benefits. Under the Equivalence Principle, the expected loss of an insurance policy to the insurer is assumed to be zero. Therefore, the expected present value of future premium (PVFP) income paid by the insured is equal to the expected present value of future benefit (PVFB) and the expected present value of future expenses (PVFE) paid by the insurer, as shown as follows: The expression of the expected PVFB, expected PVFP and the expected PVFE is given by where k represents the numbers of annual premiums payable at the beginning of each year. We assume the annual expenses charges amount to 20% of the annual premiums collected. It is noted that r b represents the risk-adjusted rate of return derived from the Capital Asset Pricing Model (CAPM). It is calculated using the following expression: HR = 1 ∶ no effect HR < 1: negative associated HR > 1: positive associated.

E(PVFP) = E(PVFB) + E(PVFE).
where denotes the correlation between unemployment rate and market rate of return, which is expressed by = cov(r u ,r m )

var(r m)
. The definition for the parameters in the Equivalence Principle framework is given in detail in Table 3. The computation of the parameters is discussed in the next section.

Results and discussion
The fitting of unemployment distribution is discussed first followed by unemployment duration, the key variables influencing unemployment, and unemployment insurance premium pricing.

Distribution fitting of unemployment rate
The fitting of the probability density function is a popular approach to explaining data patterns. The histogram in Fig. 2 shows that the unemployment data fit well with positively skewed distributions. The mathematical formula of Burr distribution and several other continuous distribution functions such as gamma, lognormal, logistic, log-logistic, and Weibull distributions are shown in Table 4. The log-likelihood function is differentiated and solved numerically to obtain the parameter estimation. Log-likelihood, Akaike Information Criteria (AIC), and Bayesian Information Criteria (BIC) are used in the model selection. The results in Table 5 show that the Burr distribution has the best fit and its hazard function is then used to estimate unemployment duration. Figure 2 shows the distribution fittings for the unemployment rate and the histogram indicates the unemployment count is left-skewed. Youth unemployment (17-24 years) accounts for 60.3% of total national unemployment. It affirms the need for a greater policy focus on youth employment despite the country's overall low unemployment rate.

Estimation of unemployment duration
The unemployment duration is measured using the hazard function. The hazard function, defined as instantaneous death rate in this study, describes the condition of an unemployed person is employed at time t provided that the person is unemployed at time t-1. Figure 3 shows that the hazard rate increases at early young age of employment and then decreases gradually after age 20. Figure 4 shows the computed results of unemployment duration.
The estimated unemployment duration shows a declining pattern with increasing age. The decline in unemployment duration with rising age is consistent with the observation that work experience is key to returning quickly to the workplace. Another plausible reason is that a person who is retiring or close to retiring will consider taking up any offer, either on part-time or contract employment (Ministry of Human Resources, 2019). The relatively young retirement age which was raised from 55 to 60 years in 2000 compared to an average of 65 years in developed countries, in tandem with rising life expectancy, are identified as the major drivers of the labor force working past the mandatory retirement age.
Further motivating this trend is the relatively young demographic profile of the country and its steady economic expansion, two major factors underlying the strong demand for skilled and experienced staff. The shorter unemployment duration for the older labor force is attributable to a combination of supply and demand factors. The employability of experienced staff also explains in part the puzzle of country's high graduate unemployment despite its low national unemployment rate.
The age 20 cohort shows the highest hazard rate estimated at 0.1457. Annualizing it (365 days), the unemployment duration is estimated at 53 days. This means that for the age 20 cohort, the maximum duration of unemployment is about 2 months. The result highlighted the contribution of this study. The youth unemployment issue merits policymakers' attention as crime and unemployment are positively associated (Fougere, et al. 2009). The youth unemployment is not unique in Malaysia. Similar patterns are observed in many countries in the world. Ciuca and Matei (2011) found that unemployment duration of 5.13 months for the less than 25 year age group. In Korea, Lim and Lee (2019) reported that 70% of the youths take 8 months or less to find a job after their graduation. In contrast, Danacica (2015) analyzed the (re) employment of young generation in Romania. For the three age groups of 15-19, 20-24, and 25-29 years examined, the author found that the chances to exit unemployment increase for the young generation below 30 years old. The unemployment duration estimates are a novel contribution to policy discussion in Malaysia. The government currently assumes a maximum of 6-month unemployment duration for allowance payment under the employment insurance scheme (EIS) that was launched in January 2018. The estimated unemployment duration indicates that the EIS may not be reflective of the current unemployment situation.

Factors influencing unemployment
The factors influencing unemployment are explored using Cox proportional hazards model. The regression results are summarized in Table 6. The variances of the coefficients (covariance) are small for all estimates, indicating low variability from the average. The coefficient estimate is used to calculate the hazard rate represented by e . The p value indicates the significance of the effects of the estimated coefficients.
From Table 6, it is observed that all age groups are found to be statistically significant in affecting national unemployment. We consider the hazard ratio to interpret the condition for the change from being unemployed to employed. For instance, the 31-36 year age group has a 64% probability to be employed. The hazard ratio of the 25-30 years age group exceeds almost 184%. This means that an  Fig. 3 Hazard function of the unemployment rate unemployed individual in this age group can quickly get a job. For the age group less than 18 years, the likelihood of being hired is only 6%, and the unemployment is contributed by the category which has less than 6-month certificate in educational attainment, or those who took general programs in their field of study. The statistical results suggest that the youths possessing full certificate or specific study program have lower unemployment compared to those equipped with general programs in their field of study. The low probability can be attributed to a majority of this age group still in school. The 19-24 years cohort has approximately 37% probability to be employed, while the 31-36 year cohort has a lower hazard ratio that could be attributed to longer working experience. The results affirm the government's concern over the relatively high youth unemployment. The policy response includes the provision of wage incentives to employees who have been out of work for a year and hiring incentives to employers to encourage youth employment. The government has also intensified skill enhancement programs to help unemployed workers return to the workplace. The current unemployment insurance scheme that pays 6 months of the allowance paid to the unemployed according to the respective percentage is considerably longer than that found in this study.
Malaysia is not alone in facing the unemployment issue. Similar unemployment condition has been observable over the past 4 decades. Katz (1974) reported that for the age cohort of 25-34, the length of unemployment (in weeks) increases in tandem with the decline in the years of schooling. It can be noticed that for the years of schooling of less than 8, the length of unemployment is 18.6 weeks. By contrast, it is about 14.2 weeks of unemployment for those with more than 12 years of schooling. For the age cohort of more than 55 onwards, the unemployment length does not vary with the respective schooling years. Those with more than 12 years of schooling averaged 21.6 weeks of unemployed which is about the same for those with less than 8 years of schooling at 22.6 weeks. Hoffman (1991) found that the age group of 19-24 has the greatest impact on unemployment. The results of the present study are consistent with the cited studies in demonstrating the importance of the age variable in influencing unemployment patterns.
On the gender effect, the p value for male is statistically significant, indicating that the male is impacted more by unemployment. This is attributable to the lower labor force participation rate of females as a significant number opt to remain as full-time housewives or are not seeking employment due to their husbands' ability to provide for the family. Some female employees work part-time as they have to care for their children. The part-time employees are classified as outside the labor force or as unemployed. As shown in Table 1, the share of females outside the labor force is 26.9%, which is more than double compared to male at only 10.4%. Further supporting the greater impact of unemployment on males is the higher hazard rate at 0.5322 compared to females at 0.0387, respectively. The figures confirm that the male is easier to get a job compared to female by almost 50%. In Russia, married women have significantly longer unemployment duration than married men (Foley, 1997). In Turkey, the report by Tansel and Tasci (2010) shows that the duration dependence of the exit rate from unemployment is different for men and women. The study found a U-shaped duration dependence for men, while for women, no duration dependence was found. For women, the percentage of leaving unemployment is lower than for men. A similar finding is obtained for Malaysia. In Finland, the women who are keen to look for employment are related to the number of children they have (Gonzalo and Saarela 2000). An increase in the number of children corresponds to a decrease in the women's exit rate into employment. A study on Russia and Ukraine (Kupets 2006) found that children do not affect the women's unemployment duration. By contrast, a different pattern is observed in the US where female has a shorter length of unemployment and a higher possibility to exit unemployment compared to male (Hoffman 1991).
On the influence of location, the results show that regional effects on unemployment are statistically significant but less pronounced compared to the level of education. An interesting outcome is that Sabah state with the highest unemployment rate is found to have an 11% higher chance to end unemployment compared to other states in Malaysia. It can be explained by the development programs mounted by the state government to enhance employability leading to a wider involvement of undergraduates and graduates from institutions of higher learning. Another plausible factor is the state's smaller population relative to its land mass and the size of economic activity.
The influence of region on unemployment has been reported in several other studies. In Romania, Danacica and Paliu-Popa (2017) found that the West Region has the shortest median survival time, while South-Oltenia Region has twice the median survival time. The study also found that the median survival duration is 754 days for rural areas. In urban, the median survival duration takes up 443 days. In Romania, Cuica and Matei (2011) show that the average duration of unemployment by counties is 6 months.
Due to its multiethnic composition, ethnicity in Malaysia is an important variable in explaining the variation in the unemployment rate and duration. The percentage to exit unemployment for the majority Malay ethnic group is estimated at 8.9% which is substantially higher than the ethnic Chinese (3.4%) and ethnic Indians (0.01%). There are limited studies on the race effect. Hoffman (1991) investigated the unemployment exit rate of black and white Americans. He found that white females have the highest chance of leaving unemployment.
Besides age, gender, region, and ethnicity, existing research has identified education attainment as a key explanatory factor for unemployment, although the results vary depending on the prevailing economic conditions. In developed countries such as the United Kingdom, Belgium, and Ireland, higher education levels shorten the unemployment duration. The chance to leave unemployment for those with tertiary education is estimated at approximately 72%. However, in East European countries as found in Romania, the higher the education level, the more difficult it is to find a job. In Greece and Spain, similar studies have reported that educational level has fewer effects on the expected unemployment duration, attesting to the influence of other variables such as economic structure and conditions in influencing employability.
In this study, education attainment across all levels is found to have a statistically significant effect on unemployment. Those with no formal education (no certificate) have a greater impact compared to those who received formal education. Among those with formal education, the lower secondary certificate category is more unlikely to exit unemployment. Interestingly, a diploma certificate stands a better chance to leave unemployment compared to a degree level, reflecting the strong industry need for technical and vocational trained rather than university graduates. On the influence of the field of study, it is observed that the labor force with agriculture and veterinary training has a high 80% probability to find employment, followed by those in education and services, respectively. The fields of social sciences, business, and law have about half the possibility to leave unemployment, rendering them with the lowest hazard ratios in the study. In Europe, Núñez and Livanos  2010) found that the disciplines that are most effective in exiting unemployment are health and welfare, education, and engineering.

Model-based unemployment insurance premium pricing
This analysis is aimed at motivating enhancements to the country's recently launched national unemployment insurance scheme. The statistical estimates derived earlier are used to formulate a model-based unemployment insurance premium. The Equivalence Principle has been used to compute the premium pricing. See Table 7 for the estimated parameters in the Equivalence Principle calculation. The current EIS has been implemented since the year 2018. Only employees earning not more than RM4000 a month are eligible for the mandatory coverage. Premium payments amounting to 0.2% of the monthly basic salary are collected from employees and employers. As shown below, this coverage can be extended based on the statistical parameters derived from the unemployment analysis. The proposed premium schedule is calculated based on the premium pricing framework described earlier. We calculate the risk-free rate of interest and market rate of return based on the median. For the unemployment rate, the first moment (expected mean) of Burr distribution is used. Table 8 shows the monthly allowance payment during the unemployment duration. We proposed an unemployment duration of a maximum of 3 months based on our findings. Some assumptions underlying the proposal include: 20% of the premium is charged for the expenses The chargeable monthly salary is constant 3 unemployment spells are allowed for the insured CAPM is considered for the expected risk-adjusted return on premiums The mortality rate is not considered during the period of insurance coverage. Figure 5 shows the comparison between the proposed premium and the current EIS. Based on the current EIS, it is seen that a worker with a monthly salary of RM1,000, the worker pays an annual premium of RM24.00. The computed schedule shows an annual premium of RM7.95. The wide pricing gap highlights the need to examine the underlying cost structure including profit margin, reserves accumulation, and claims-paying ability. The salary proportion is now reduced to 0.066% compared to 0.2% and there is no cap on eligibility based on wage level. The median household income has been taken into consideration in the proposed scheme. The premium paid by the respective household group is also lower compared to the existing scheme. The bottom 20% income group (B20) pays RM23.85 for the annual premium for an average monthly household income of RM3000. The annual premiums for middle-income (M40) and top income (T20) groups are computed at RM49.30 and above RM80.00 per annum, respectively. Figure 6 displays the household group income in Malaysia as of the year 2019. The detailed premium schedule is given in Supplementary Table 15.

Study implications
The findings highlighted that age has the most influence on unemployment in Malaysia among the covariates investigated in this study. The youth group (19-30 years) is found to be the most vulnerable given its higher unemployment rate and longer unemployment duration. To address this key unemployment issue, the government has launched various initiatives such as job and skills matching and entrepreneurship, including seeking employment and entrepreneurial opportunities in business start-ups and small-and medium-sized enterprises. In the 2022 budget announcement, the Ministry of Finance Malaysia allocated RM4.8 billion to fund the job guarantee program (JaminKerja) to create 600,000 job opportunities. The budget will also target 220,000 trainees for various training and upskilling programs with an additional allocation of RM1.1 billion. On reskilling and upskilling, a dedicated statutory body, the Human Resource Development Corporation (HRD Corp) has launched a digital platform called e-LATiH that will provide them access to more than 300 online courses to meet the demands of various industries. In addition, e-LATiH also provides a job matching platform to job seekers and employers. Given the crucial role of the young generation in nation-building, these initiatives are important to address the social ills and underutilization of human capital due to high youth unemployment as affirmed in this study.
The current unemployment insurance scheme implemented in Malaysia considers charges of 0.2% from employers and 0.2% from employees respectively, with the coverage limited to employees whose monthly salary is not more than RM4000 a month. This is a major drawback of the current EIS. Based on the statistical inference, the unemployment scheme derived in this study requires a lower insurance premium and no imposition of salary limits, thereby allowing for universal coverage of the country's entire workforce. The proposed insurance scheme is, therefore, more reasonable that can be extended to every household in the country.

Conclusion
Unemployment is an overarching development issue regardless of a country's income level and a crucial policy challenge irrespective of the prevailing level and duration of joblessness. This study is the first extensive analysis using statistical tools to investigate the determinants of unemployment in Malaysia, an upper-middle-income country. It models the distribution of unemployment and provides statistical support for extending the coverage of the existing unemployment insurance scheme. The Burr distribution is used to fit the unemployment data and the Cox proportional hazards model is applied to determine the key influencing variables. The study sheds new light on the challenges facing Malaysia despite its low unemployment situation. Based on the statistical inference, a feasible unemployment insurance scheme is developed to provide wider coverage at a lower cost to employees. The study found that the 19-24 year age group has the highest unemployment rate, affirming the need for the youth unemployment issue to take center stage in policy formulation as well as in enhancing the design of the existing national unemployment insurance scheme. The statistical evidence shows that it is easier for the male cohort to exit unemployment than for women, suggesting that greater policy focus on female unemployment will have a greater overall impact on development. The study also found significant differences across states, whereby Sabah has the highest rate of unemployment, but further analysis shows that the regional effect is less pronounced than education level. Consistent with the findings in other countries, those with higher education have a better chance to exit unemployment. Based on fundamental design principles, the unemployment insurance premium framework applied in the study suggests that the existing scheme can be expanded with lower premiums charged and without any cap on wages for employees to be eligible for coverage. The proposed insurance scheme framework developed in this paper can be used to motivate enhancements to the existing scheme after considering other factors such as age, education level, and field of study in a dynamic risk-based insurance scheme. To design a more flexible risk-based unemployment insurance scheme, the significant covariates identified in the Cox regression hazard model can be used as a guide to incorporate other risk factors.

Consent for publication
The data are provided by the national statistical agency, Department of Statistics Malaysia (DOSM). Citation has been done to acknowledge the contribution.


s2fieldsofstudy: Economics Mathematics

text:
Uncovering a two-phase dynamics from a dollar exchange model with bank and debt

We investigate the unbiased model for money exchanges with collective debt limit: agents give at random time a dollar to one another as long as they have at least one dollar or they can borrow a dollar from a central bank if the bank is not empty. Surprisingly, this dynamics eventually leads to an asymmetric Laplace distribution of wealth (conjectured in [22] and shown formally in a recent work [18]). In this manuscript, we carry out a formal mean-field limit as the number of agents goes to infinity where we uncover a two-phase (ODE) dynamics. Convergence towards the unique equilibrium (two-sided geometric) distribution in the large time limit is also shown and the role played by the bank and debt (in terms of Gini index or wealth inequality) will be explored numerically as well.

Introduction
Econophysics is a subfield of statistical physics that apply concepts and techniques from traditional physics to economics or finance [12]. One of the primary goal of this area of research is to explain how various economical phenomena could be derived from universal laws in statistical physics under certain model assumptions, and we refer to [15] for a general review.
There several motivations for the study of econophysics models: from the perspective of a policy maker, how to influence/control the emerging wealth inequality (measured by Gini index) in order to mitigate the alarming gap between rich and poor is a central issue to be dealt with. From a mathematical viewpoint, the fundamental mechanisms behind the formation of macroscopic phenomena, for instance various possible wealth distributions from different agent-based money exchange models, must be thoroughly understood. For a given (stochastic) agent-based model, we would like to identify a limit (deterministic) dynamics when we send the number of individuals to infinity, and then the deterministic system will be further analyzed with the intention of proving its convergence to equilibrium (if there is one) for large time. This paradigm has been implemented in vast amount of works across different fields of applied mathematics, see for instance [1,6,7,20].
In this work, we consider a simple mechanism for money exchange involving a bank, meaning that there are a fixed number of agents (denoted by N ) and one bank. We denote by S i (t) the amount of dollars the agent i has at time t and we suppose that N i=1 S i (0) = N µ for some fixed µ ∈ R + so that N µ ∈ N + . Moreover, we denote by B * := B c + B d the initial amount of dollars in the bank, where B c and B d represent the amount of dollars owned by the bank in the form of "cash" and in the form of "debt" (borrowed by agents), respectively. Also, we introduce another parameter ν ∈ R + so that N µ ν ∈ N + and set B * = N µ ν.
The model investigated in this work was proposed in [22] and revisited by [18]: at random time (generated by an exponential law), an agent i (the "giver") and an agent j (the "receiver") are picked uniformly at random. If the "giver" i has at least one dollar (i.e. S i ≥ 1) or if the central bank has "cash" (i.e. B c ≥ 1), then the receiver j receives a dollar. Otherwise, when the receiver i has no dollar and the bank has no cash, then nothing happens. We illustrate the dynamics in figure 1 explaining the three cases when agent i is picked to give one dollar to agent j. From now on, we will call this model the unbiased exchange model with collective debt limit and it can be represented by (1.1).
Notice that when the bank gives a dollar to agent j, there is still one dollar withdrew from the giver i, i.e. the debt of agent i increases (represented in red in figure  1). The debt of agent i could be reduced once it will become a "receiver". It is also important to notice that in this model the bank never losses money, it just transforms its "cash" B c into "debt". Without the bank, agents can only give a dollar when they have at least one dollar (hence no debt is allowed). In this case, the model is termed as the one-coin model in [16], the unbiased exchange model in [2,5], and the mean-field zero range process in [19].

Remark.
In order to have the correct asymptotic as the number of agents goes to infinity N → +∞, we need to adjust the rate λ by normalizing by N so that the rate of a typical agent giving a dollar per unit time is of order 1.
The fundamental question of interest is the exploration of the limiting money distribution among the agents as the total number of agents N and time t become large. To foresee the behavior of the dynamics under these limits, we perform a numerical simulation with N = 10, 000 agents over 10, 000, 000 steps. In figure  2-left, we plot the evolution of the cash in the bank B c over time. We observe a first phase where the cash B c is decaying linearly until it reaches zero (after roughly 2 millions steps). After that, there is a second phase where B c remains close to zero. In figure 2-right, the distribution of the wealth distribution is plotted after 10 millions steps. This distribution is well approximated by an asymmetric Laplace distribution given in (1.4).
This numerical result will be explained by our derivation in section 2 and the following asymptotic analysis in section 3. Our approach illustrated by figure 3 consists in two steps. Our first step consists in deriving the limit dynamics as the number of agents goes to infinity N → +∞ (section 2). With this aim, we introduce the probability distribution of wealth: with p n (t) = {"probability that a typical agent has n dollars at time t"}. The evolution of p(t) will be given by a (deterministic) nonlinear ordinary differential equations. To fully justify this transition from a stochastic interacting systems into a deterministic set of ODEs, one needs the so-called propagation of chaos [21]. We do not investigate the proof in this manuscript, the derivation has been rigorously justified in various models arising from econophysics, see for instance [2][3][4][5]10,11,13].  Figure 1: Illustration of the unbiased exchange model with collective debt limit: at random time, a "giver" i has to give one dollar to a "receiver" j. Three cases are possible. Case 1: if the "giver" i has at least a dollar (i.e. S i ≥ 1), then it gives it to j. Case 2: if i does not have a dollar (S i ≤ 0) and the central bank has cash (i.e. B c ≥ 1), then j receives one dollar from the bank and the debt of i is increased by one. Case 3: the giver i and the bank do not have any money, nothing happens.  The additional difficulty here is that the evolution of p(t) is split into two phases. Indeed, the evolution changes at the first time when there is no more cash in the bank. We will denote by t * the time at which such event occur, i.e. at t * the bank is empty for the first time. The evolution equation of p(t) takes the form of: where the exact expressions for the operator Q 1 and Q 2 are given by (2.5) and (2.15), respectively. Our second step is to investigate the asymptotic behavior of the probability mass function p(t) as t → +∞, more precisely its convergence toward an asymmetric Laplace distribution given in (1.4): where the three positive parameters ρ 0 , α and β only depend on µ (average wealth per agent) and ν (the ratio between the bank and agent's combined wealth). Our proof relies on the so-called entropy method which gives a rigorous proof of the convergence of p(t) toward the asymmetric Laplace distribution ρ, see section 3.2. Moreover, a standard linearization analysis performed in section 3.3 allows us to obtain an exponential decay result for the linearized entropy (near equilibrium), provided that the parameters µ and ν fulfill certain criteria. Finally, we would like explore the role played by the bank in terms of wealth inequality. In particular, we conjecture that the inclusion of the bank and possibility for agents to go into debt will lead to accentuation of wealth inequality (measured by the Gini index), compared to the usual unbiased exchange model (without the presence of a bank). This type of Gini index comparison conjecture (shown numerically in section 4) makes sense at least heuristically but we are not able to provide a rigorous proof.
Although we will only investigate a specific binary exchange models in the present work, other exchange rules can also be imposed and studied, leading to different models. To name a few, the so-called immediate exchange model introduced in [14] assumes that pairs of agents are randomly and uniformly picked at each random time, and each of the agents transfer a random fraction of its money to the other agents, where these fractions are independent and uniformly distributed in [0, 1]. The so-called uniform reshuffling model investigated in [3,12,17] suggests that the total amount of money of two randomly and uniformly picked agents possess before interaction is uniformly redistributed among the two agents after interaction. For closely related variants of the unbiased exchange model we refer to the recent work [2]. For other models arising from econophysics, see [4,8,9,18].

Formal mean-field limit
In this section, we carry out a formal mean-field argument of the stochastic agentbased dynamics (1.1) as the number of agents N goes to infinity. Even our analysis is not completely rigorous, the resulting system of ODEs (i.e., (2.4)-(2.5) in Phase I and (2.14)-(2.15) in Phase II) admits a unique equilibrium distribution as t → ∞, which is a two-sided geometric distribution and can be well-approximated by an asymmetric Laplace distribution when µ 1 (i.e., when the initial average amount of dollars per agent becomes large). As our predication of the limiting distribution of money based on the large time behavior of the ODE system (2.14)-(2.15) coincides with the conjecture in [22] and the work in [18], it strongly indicates that our heuristic mean-field analysis actually captures the correct deterministic behavior of the underlying stochastic dynamics when we send N → ∞.
The unbiased exchange model with collective debt limit can be written in terms of a system of stochastic differential equations, thanks to the framework set up in [2] for the study of the basic unbiased exchange model as well as some of its variants. Introducing {N (i,j) t } 1≤i,j≤N , which represents a collection of independent Poisson processes with constant intensity λ N , the evolution of each S i is given by: t "i gives to j" where we use the notation δ 0 (B c ) := 1{B c = 0}. By the obvious symmetry, we can focus on the case when i = 1 and notice that whenever B c ≥ 1, the SDE for S 1 simplifies to If we introduce then the two Poisson processes N 1 t and M 1 t are of intensity λ := 1. Motivated by (2.2), we give the following definition of the limiting dynamics of S 1 (t) as N → ∞ from the process point of view, providing that B c ≥ 1.

Definition 1 (Mean-field equation -Phase I)
We define S 1 to be the compound Poisson process satisfying the following SDE: in which N 1 t and M 1 t are independent Poisson processes with unit intensity. We denote by p(t) = (. . . , p −n (t), . . . , p −1 (t), p 0 (t), p 1 (t), . . . , p n (t), . . .) the law of the process S 1 (t), i.e. p n (t) = P S 1 (t) = n . Its time evolution is given by the following dynamics: d dt

Remark.
It is readily seen that p(t) is exactly the law of a continuous-time symmetric simple random walk on Z defined via 1} being a sequence of independent Rademacher random signs and N t being a Poisson clock running at unit rate.
The distribution p(t) naturally preserves its mass (since it is a probability mass function) and the mean value (as the total amount of money in the whole system is preserved), thus if we introduce the affine subspaces: and S + µ = {p ∈ S µ | p n = 0 for n < 0}, then it is clear that the unique solution p(t) of (2.4)-(2.5) with p(0) ∈ S + µ satisfies p(t) ∈ S µ for all t > 0. Moreover, if we define the average amount of "debt" per agent as this quantity will be non-decreasing: d dt D[p(t)] = p 0 (t) ≥ 0. Since the average amount of debt each agent can sustain in the underlying stochastic N -agents system is at most µ ν, we therefore terminate the evolution of Phase I (2.4)-(2.5) until

Remark.
Thanks to the previous remark, E|X t | t→∞ − −− → +∞. By the identity |X t | = X + t + X − t and the obvious symmetry, we deduce that is also unbounded as t → ∞. This observation ensures the finiteness of the t * .
At the level of the agent-based system, after the first time when there is no cash in the bank, i.e., when t ≥ t stoc * := min{τ > 0 | B c (τ ) = 0}, the analysis is much more involved so heuristic reasoning plays a major role in this manuscript. We notice that we have the following basic relations for all time t ≥ 0: and B c ≥ 0. Therefore, the evolution of B c is much faster than the evolution of each of the S i 's, indicating that (2.10) is really a "fast" dynamics compared to (2.1). These observations motivate the next definition:

Definition 2 (Mean-field equation -Phase II)
We define S 1 (for t ≥ t stoc * ) to be the compound Poisson process satisfying the following SDE: 11) in which N 1 t and M 1 t are independent Poisson processes running at the unit intensity. Moreover, Y = Y (t) ∼ B (r(t)) and Z = Z(t) ∼ B (q 0 (t)) are independent Bernoulli random variables with parameters r := P S 1 ≥ 1 (2.12) and respectively.
Remark. We will illustrate the basic intuition behind Definition 2. Under the large population limit N → ∞, we denote by It is easily seen that the transition rates of B c can be described by figure 4 below.
(rate ) (rate ) That is, the transition B c → B c + 1 occurs when a "rich" agent is being picked (S i ≥ 1) and gives a dollar to an agent in "debt" (S j < 0). Similarly, the transition B c → B c − 1 happens when an agent without dollar is being picked (S i ≤ 0) and gives a dollar to an agent with debt (S j ≥ 0) and moreover B c ≥ 1. As the evolution of B c is a "fast" dynamics, we assume its distribution will relax to its ergodic invariant distribution within a time-scale that is negligible compared to the evolution of each of the S i 's. Thus, from the following detailed balance equation at equilibrium q n · r · d = q n+1 · (r + p 0 ) · (d + p 0 ), one arrives at q n = r d (r+p 0 ) (d+p 0 ) n q 0 for all n ≥ 0. Taking into account that n≥0 q n = 1, we end up with "in debt" "rich" (rate ) when (rate ) when which is exactly (2.13). On the other hand, the transition rates of S 1 can be described by figure 5 below.

Large time behavior
Although our mean-field analysis is not completely rigorous, we will soon see that our system of ODEs (i.e., (2.4)-(2.5) for t ≤ t * and (2.14)-(2.15) for t ≥ t * ) converges to a double-geometric distribution as t → ∞, which resembles an asymmetric Laplace distribution when µ becomes large. Taking into account the conjecture in [22] and the work of [18], as well as several numerical experiments, we believe that the two phase dynamics (1.3) accurately captures the mean-field behavior of the stochastic N -agents dynamics as N → ∞.

Elementary properties of the limit equation
In Phase II, the collective debt d (2.16) is preserved, therefore we can further restrict the affine space where the solution p(t) of (2.14)-(2.15) lives: it is straightforward to check that if p(t * ) ∈ S µ,ν then p(t) ∈ S µ,ν for all t ≥ t * . Now, we identify the unique equilibrium solution associated with (2.14)-(2.15) in this space.

Proposition 3.1
The unique equilibrium solution of (2.14)-(2.15) in S µ,ν , denoted by p * = {p * n } n∈Z , is given by in which Proof. From the evolution equation in Phase II (2.14)-(2.15), it is straightforward to check that the equilibrium distribution takes the ascertained form (3.2). As we also have p * ∈ S µ,ν , we impose that n∈Z n p * n = µ, and − n≤0 n p * n = µ ν.
These constraints lead us to Combining these relations with the elementary identity r * + p * 0 + d * = 1, we obtain the desired result.
Remark. Under the settings of Proposition 3.1, if we take µ 1 while keeping ν fixed, the two-sided geometric distribution p * (3.2) can be well-approximated by the asymmetric Laplace distribution given by (1.4) and conjectured in [18,22]. Furthermore, we have an estimation of the parameters ρ 0 , α, β of this distribution with Indeed, by taking µ 1 in (3.3), we can perform the following simple asymptotic analysis for p * 0 : Similarly, the asymptotic behavior of α can be seem from the asymptotic analysis for log 1 r * r * +p * 0 . Thanks to (3.3), (3.4) and the asymptotic argument for p * 0 above, we have The argument for β is pretty similar and hence will be omitted here.

Convergence to asymmetric Laplace distribution
The main goal of this section is to prove the convergence of the solution p(t) in Phase II (2.14)-(2.15) to its two-sided geometric equilibrium solution p * as t → ∞. The essence of the method consists in studying the (relative) entropy dissipation between the two distributions. We remind the formula for estimating the relative entropy: Theorem 1 Let p(t) = {p n (t)} t≥t * be the unique solution to (2.14)-(2.15) with p(t * ) ∈ S µ,ν , then for all t ≥ t * we have Consequently, if n≥0 p n (t * ) c n + n≤0 p n (t * ) c −n < ∞ for some c > 1, then the solution p(t) converges strongly in p for 1 < p < ∞ as t → ∞ to the two-sided geometric distribution p * .

Proof.
We split the computation of p n+1 log p * 0 + n log and On the other hand, straightforward computations lead us to Combining (3.11) and (3.12) gives rise to the advertised entropy dissipation result (3.7). The last statement of Theorem 1 can be handled in a pretty similar way as in the case for the basic unbiased exchange model (without the presence of a bank), see for instance [2] or [19].
To illustrate the convergence of p(t) toward the equilibrium p * , we run a simulation and plot the evolution of p(t) at different times (see figure 6-left) as well as the evolution of the relative entropy D KL (p(t) || p * ) over time (see figure 6-right). Notice that the decay of the relative entropy changes abruptly around t ≈ 200 which corresponds to the transition of the dynamics from phase I to phase II. During the phase II, the decay of the relative entropy is well-approximated by a decay of the form: (3.13) where the coefficients c 1 = .674 and c 2 = .182 are estimated using mean-square error criteria. We emphasize here that the ansatz for the decay of the relative entropy D KL (p(t) || p * ) (3.13) is inspired from a very recent work on the vanilla unbiased exchange model [5], which corresponds to the special case of the model at hand where ν = 0.

Linearization analysis
We perform a standard linearization analysis near the two-sided geometric equilibrium distribution p * . The hope is that for some ranges of parameter choices for µ and ν, the linearized entropy will decay exponentially fast in time. As we will soon wealth ($) probability equilibrium time (t) fitting: relative entropy relative entropy realize, to enlarge the parameter choices of µ and ν for which the linearized entropy decreases exponentially fast, one has to find the best possible constants for certain Poincaré-type inequalities. We linearize the ODE system (2.14)-(2.15) near its equilibrium p * , i.e., we set p n = p * n + ε w n for 0 < ε 1 for all n ∈ Z to obtain the following linearized system around p * : where we set It is straightforward to check that for all t. Moreover, if we denote the linearized entropy associated with the linearized system (3.14) by a direct computation yields We can bound α 2 as follows. Introduce {α n } n =0 by defining α n = 1 r * for n ≥ 1 and α n = 1 r * +p * 0 for n ≤ −1. Similarly, let {λ n } n =0 be defined by λ n = γ 1 for n ≥ 1 and λ n = γ 2 for n ≤ −1, where γ 1 and γ 2 are constants whose choices will be optimized later. The definition of α, together with the fact that w ∈ V, allows us to deduce that Optimizing the values of γ 1 and γ 2 gives rise to

Thus, we obtain
In a similar fashion, we also get as well as Thus, we have the following upper bound: (3.20) On the other hand, we have the following simple lower bound on (3.21) To alleviate the writing, we define . (3.28) Then we arrive at the following differential inequality for the linearized entropy E[w]: If we recall the upper bound (3.19) on , in order to have an exponential decay result for the linearized entropy, it suffices to ensure that For convenience, we denote by G µ,ν the collection of the parameter values (µ, ν) ∈ R 2 + for which the condition (3.30) is satisfied. We claim that this set is non-empty, whence the linearized entropy will decay exponentially fast in time for all (µ, ν) ∈ G µ,ν .

Proof.
To demonstrate the non-emptiness of the set G µ,ν it suffices to present a specific example. For this purpose, we can take (µ, ν) = (0.01, 0.001) to obtain The second part of the proposition follows immediately from the definition of G µ,ν and the differential inequality (3.29).

Remark.
The lower bounds (3.21) and (3.22) are far from being optimal as their arguments rely sole on some elementary Cauchy-Schwarz inequalities and the constraint that w ∈ V has not been exploited at all. We speculate that these lower bounds (Poincaré-type inequalities) can somehow be significantly sharpened, leading to a refined version of Theorem 2.

Debt induced wealth inequality
We would like to conclude our analysis of the model by investigating the effect of the bank on the inequality of the wealth distribution, i.e. does the bank increase or decrease inequality? To measure inequality, we use the Gini index G which is usually an indicator between 0 (no inequality) and 1 (maximum inequality) and measures the inequality in the wealth distribution, and a higher Gini index indicates worse equality in the society (in terms of the distribution of wealth).

Definition 3 (Gini index)
For a given probability mass function p ∈ P(Z) with mean µ ∈ R + , the Gini index of p is defined by Equivalently, a probabilistic definition of G[p] is given by where S has p as its probability mass function and S is an independent copy of S.
Remark. Without debt, i.e. if p n = 0 for n < 0 or S ≥ 0, one can show that the Gini index G is always between 0 (no inequality) and 1 (maximum inequality). Indeed, using triangular inequality: However, the Gini index is no longer bounded if the wealth distribution includes debts. For instance, taking a Gaussian distribution with mean µ > 0 and variance σ 2 , i.e. X ∼ N (µ, σ 2 ), we obtain: Thus, if σ µ, the Gini index could be arbitrary large.
We would like to compare the evolution of Gini index for the basic unbiased exchange model and the model at hand (with bank and debt) at the level of the deterministic ODE system, starting from the same initial distribution. First, we recall that the unbiased exchange model investigated in [2,13,16,19] is a special case of the model studied in this work with ν = 0, meaning that the bank does not exist and agents are not allowed to go into debt. The rigorous mean-field analysis shown in [2,13,19] implies that, if we denote by q(t) = (q 0 (t), q 1 (t), . . .) the law of the amount of dollars a typical agent has as N → ∞, its time evolution will be given by: with Q[q] n := q n+1 +r q n−1 − (1 +r) q n , for n ≥ 1, q 1 −r q 0 for n = 0, (4.7) wherer := n≥1 q n . Intuitively, the model with bank investigated in this manuscript permits agents without dollars in their pocket or agents in debt to be picked to give, thus we speculate that the Gini index of the distribution p(t) solution of (1.3) is always larger than the corresponding Gini index of the distribution q (if they start from the same initial condition).
We provide in figure 7-left the evolution of the Gini index for different values of ν (which measures the ratio of the initial total amount of cash in the bank to the agents' combined initial wealth). Without the bank, i.e., when ν = 0, the Gini index quickly reaches its maximum around G ≈ .8. However, when the bank comes into play (ν ≥ .25), the Gini index could eventually exceed 1. Moreover, the phases I and II are clearly seen in the evolution of the Gini index. In phase I, the Gini index grows like O( √ t) which is due to the diffusive nature of the dynamics coupled with formula (4.5). In phase II, the growth of the Gini index starts to saturate and the curve converges to its maximum. In figure 7-right, we provide the corresponding wealth distribution at the final time of the computation, i.e. p(t = 5000), depending on the wealth of the bank ν. Without bank (ν = 0), the distribution is exponential but it becomes an asymmetric two-sided exponential for ν > 0. Notice that in the case where ν = 5 (which means that the bank has 5 times more resources than the agents), the dynamics has not reached yet phase II at time t = 5000 as we observe in the evolution of the Gini index G in figure 7-left. As a result, the distribution p(t = 5000) for ν = 5 is still far away from equilibrium and the characteristic cusp at zero dollar has not yet appeared.
Although we are not able to provide a complete proof of this natural and compelling conjecture based on heuristic reasoning and numerical simulations. The wealth ($)  that p(t) is the solution of (2.4)-(2.5) for t ∈ [0, t * ] with p(0) ∈ S + µ , then the Gini index of p(t) is non-decreasing with respect to time. Proof. We resort to the probabilistic representation of the Gini index (4.2), so we introduce independent random variables S and S having p as their common probability mass function. If we set Z = |S − S | and let z = {z n } n≥0 to be law of Z, it is trivial to see that the evolution of z satisfies d dt z n = 2 (z n+1 + z n−1 − 2 z n ) for n ≥ 1, 2 (z 1 − 2 z 0 ) for n = 0, (4.8) from which we deduce that d dt In other words, the phase I evolution tends to accentuate the wealth inequality and the proof of Proposition 4.1 is completed.

Conclusion
In this manuscript, the so-called unbiased exchange model with collective debt limit is investigated, which can be viewed as an extension of the vanilla unbiased exchange model proposed in [12] and revisited in [2,16,19]. Although the inclusion of bank creates additional difficulty in the analysis and only a formal mean-field argument is presented, we found that the prediction of the asymptotic distribution of wealth as N → ∞ and t → ∞ based on our two-phase dynamics (1.3) agrees with the results reported/conjectured in earlier work [18,22]. To the best our of knowledge, there are no (even formal) mean-field analysis for the model at hand prior to the present work and we believe that our work leaves many open and challenging questions to be investigated on a more rigorous ground. For instance, is it possible to provide a rigorous proof of the propagation of chaos as N → ∞ to arrive at the two-phase dynamics (1.3) ? How can we extend the recent work [4] to justify the ansatz for the decay of the relative entropy D KL (p(t) || p * ) conjectured in (3.13) (for a fairly generic initial datum) ? Is it possible to show rigorously the monotonicity of the Gini index of p(t) (solution of (1.3)) with respect to ν during phase II ? The answer to these questions will enable us to have a better understanding of the role played by the bank and debt. Lastly, we emphasize that econophysics models involving bank and debt have not received enough attention in the past few years, and the model investigated in this manuscript can be viewed as the representative of a promising and intriguing direction for further research in the area of econophysics. For instance, it would be interesting to introduce a bank in other closely related variants of the basic unbiased exchange model, such as the so-called poor-biased or rich-biased exchange model proposed in [2].


Spatial frequency of unstable eigenfunction of the core-periphery model incorporating differentiated agriculture with transport cost

We extend the core-periphery model with differentiated agriculture with transport cost to those on discrete and continuous multi-regional space and investigate the stability of homogeneous stationary solution, especially for the model on a continuous periodic space. Unstable eigenfunctions can be stable when manufacturing transport cost is sufficiently high, but become unstable when the cost is lower than a certain critical point. As the transport cost decreases further below another critical point, the eigenfunction becomes stable again. It can be observed numerically that the unstable area, which is the range of the transport cost between those two critical points, generally tends to expand with the frequency increases limited to even or odd numbers only.

Introduction
The core-periphery model proposed by Krugman (1991) provides a basic framework for modeling economic agglomeration based on the general equilibrium theory. In its most basic form, symmetric equilibrium, where all regions are uniformly populated, is stable when transport cost is sufficiently high, while it becomes unstable when the transport cost is sufficiently low, producing a spatially inhomogeneous structure.
In the most basic formulation of the model, agricultural goods are assumed to be undifferentiated and, in addition, their transport is assumed to have no cost. Meanwhile, Fujita et al. (2001, Chapter 7) introduce a model incorporating differentiated agriculture and its transport cost 1 . This model has an interesting property that symmetric equilibrium stabilizes not only when the manufacturing transport cost is sufficiently high but also when the manufacturing transport cost is sufficiently low 2 .
We extend the core-periphery model with differentiated agriculture and its transport cost to finite or infinite multi-regional spaces and consider the stability of stationary solution, especially in the case of one-dimensional periodic continuous space. In the continuous space model, disturbances to the stationary state are decomposed into eigenfunctions, each with its own spatial frequency. The value of the spatial frequency corresponds to the number of cities or economic centers formed when workers and capitals begin to agglomerate.
Depending on the value of various parameters, there would be eigenfunctions which are stable over the entire range of manufacturing transport cost (stable eigenfunctions). On the other hand, there are eigenfunctions that become unstable at certain manufacturing transport cost values (unstable eigenfunctions). We see that unstable eigenfunctions can be stable when manufacturing transport cost is sufficiently high, but become unstable when the cost is lower than a certain threshold value of the transport cost (upper critical point). As the transport cost decreases further below another threshold (lower critical point), the eigenfunction becomes stable again. The behavior of these unstable eigenfunctions is similar to the solution of a linear problem around a homogeneous stationary solution of the original two-regional model by Fujita et al. (2001, Chapter 7).
We obtain expression for the eigenvalue corresponding to the eigenfunction of each frequency and investigate its frequency-dependent behavior. From the numerical observations, it is suggested that as the spatial frequency number of the eigenfunction increases within the even-only or oddonly range, the upper critical point would be non-decreasing (generally increasing) and the lower critical point would be non-increasing (generally decreasing), resulting in expansion of the unstable area where the eigenvalue takes positive value. Meanwhile, it is often observed that the value of lower critical point become larger with increasing frequency that are not limited to even or odd numbers. Nevertheless, if the frequency increases sufficiently, the unstable area generally expands.
1 Picard and Zeng (2005) consider a quasi-linear model that takes into account differentiated agriculture and its transport cost. In addition, Zeng and Takatsuka (2016) and Zeng (2021) even treat the case of a footloose entrepreneur model. In the context of international trade theory, models including agriculture have also been devised. For example, Takatsuka and Zeng (2012) consider the transport cost of the homogeneous agricultural sector, and Zeng and Kikuchi (2009) deal with differentiated agriculture with its transport cost.
2 Therefore, the symmetric equilibrium becomes unstable when the transport cost is moderate.
The rest of the paper is organized as follows. Section 2 derives the model. Section 3 discusses the homogeneous stationary solution and the eigenvalue of the linearized problem around the stationary solution. Section 4 examines how the eigenvalues depend on various parameters, mainly frequency. Section 5 discuses the results and gives conclusions. Section 6 gives proofs of several propositions omitted in the main part.

The model
This section provides an overview of Dixit-Stiglitz framework (Dixit and Stiglitz (1977)) in the context of the differentiated agriculture with transport cost and applies it to derive the multi-regional version of the model which we handle.

Dixit-Stiglitz framework
The economy is supposed to have a continuously infinite variety of manufactured goods and agricultural goods. Let m(i) ≥ 0 and a(j) ≥ 0 denote demands for the variety i ∈ [0, V M ] of the manufactured goods and the variety j ∈ [0, V A ] of the agricultural goods. Price of i-th variety of the manufactured good is denoted by p M (i) and that of j-th variety of agricultural good is denoted by p A (j). Given an income Y , and the prices p M (i), p A (j), each consumer faces the following budget constraint.
Facing (1), the consumer is supposed to maximize the following utility function. where Here, parameters η > 1 and σ > 1 stand for the elasticity of substitution between any two varieties of agricultural goods and that of manufactured goods, respectively. Following Fujita et al. (2001, Chapter 4), we solve this problem in two steps. First, the minimum cost to achieve the partial utility M and A is calculated, and then the optimal expenditure on M and A under the budget constraint is determined 3 .
We obtain a function of consumer demand for the j-th variety of the agricultural goods 4 as where Y is the income of the consumer, and G A is the price index of the agricultural goods defined by 2.2 Multi-regional modeling (Discrete space) Based on the framework described above, we model a two-sector spatial economy consisting of multiple regions. The modeling of manufacturing sector has already been described in detail in Fujita et al. (2001, Chapter 5 and 6), so we only mention it briefly. On the other hand, the agricultural sector is discussed in more detail here, since Fujita et al. (2001, Chapter 7) only mention the two-regional case. In this subsection, the economy consists of R discretely located regions indexed by r or s = 1, 2, · · · R. Following Fujita et al. (2001, Chapter 7), the total population of the mobile and immobile workers are normalized to be µ and 1 − µ, respectively. Then, the mobile and immobile population in region r is denoted by µλ r and (1 − µ)φ r , respectively. By these, it is easy to see that R r=1 λ r ≡ 1 and R r=1 φ r ≡ 1. Let w M r and w A r denote the nominal wage of the manufacturing and agricultural workers, respectively. Then, the total income of region r is given by The prices of any variety of agricultural goods produced in region r have equal value p A r . Assume that the transportation of the agricultural goods incurs the iceberg transport cost. That is, to transport one unit of agricultural goods from r to s, T A rs ≥ 1 times as much agricultural goods must be shipped. Then, the price in region s of a variety of agricultural good produced in region r is p A rs = p A r T A rs . The number of varieties of agricultural goods produced in region r denoted by n A r is assumed to be proportional to the agricultural population in region r as n A r = φ r , r = 1, 2, · · · , R.
This implies that the total number of varieties of agricultural goods in the economy is normalized to 1.
Then from (6) and (8), the price index of the agricultural goods in r is given by By assuming p A r = w A r , (∀r), we obtain the following price index equation Let us derive the equation for the nominal wage of the agricultural workers. The demand for the j-th variety of agricultural goods by an consumer with income Y is given by (5). Thus, the demand from region s for one variety of agricultural good produced in region r can be given as We obtain the sales of an agricultural producer in region r as the following by adding up T A rs times the amount in (10) over all regions.
Assume that the total supply of agricultural goods in region r is equal to the population of the agricultural workers there given by (1 − µ)φ r . Since the number of varieties of the agricultural goods produced there is n A r = φ r according to (8), the supply of one variety of the agricultural goods produced there is given by (1 − µ)φ r /φ r = 1 − µ. Thus, the condition that the supply 1 − µ and demand (11) for each variety of the agricultural goods are equal, together with the assumption that p A r = w A r , yields For the price index G M r and nominal wage w M r in the manufacturing sector, exactly the same equations are derived as in the one-sector racetrack model (See (5.4) and (5.5) of Fujita et al. (2001, p.63)). That is As mentioned in Fujita et al. (2001, Chapter 7), the real wage of the mobile workers is We adopt the following ad-hoc dynamics in which the mobile workers flow into regions with above-average real wage, while they flow out of region with below-average real wage. The average real wage is defined by R s=1 ω s λ s and the dynamics is described by where γ > 0 stands for adjustment speed. The equations (7), (9), (12), (13), (14), (15), and (16) constitute the following system of nonlinear algebraic differential equations.

Multi-regional modeling (Continuous space)
Let N be a n(≥ 1)-dimensional manifold in a continuous space. It is easy to extend the system (17) to a system on N . The mobile workers are continuously distributed on N and the density at time t ∈ [0, ∞), region x ∈ N is denoted by µλ(t, x) ≥ 0 which satisfies M λ(t, x)dx ≡ 1 for any t. The immobile workers are also distributed on N and the density at region x ∈ N is denoted by (1 − µ)φ(x) ≥ 0 which satisfies N φ(x)dx ≡ 1. Then, the system (17) naturally becomes the following continuous system 5 .
In the following, we consider in particular a circumference S of radius r > 0 as a manifold N and discuss the system (18) on S. The functions T A (x, y) ≥ 1 and T M (x, y) ≥ 1 represent the amount of the agricultural and manufacturing goods that must be shipped to transport one unit of the agricultural and manufacturing goods, respectively, from x to y. They are also restricted to the forms respectively, where τ A ≥ 0 and τ M ≥ 0 are cost parameters. Here, the distance function |x − y| is defined as the shorter distance between x and y along the circumference S. Then, the parameters τ A , τ M , η, and σ often appear in the form τ A (η − 1) and τ M (σ − 1), so it is convenient to define In the following, we identify a function on S with a corresponding periodic function on the interval [−π, π], that is, , ω(t, x), and λ(t, x) of (t, x) ∈ [0, ∞) × S are identified with the corresponding functions which are periodic with respect to spatial variables Y (t, θ), G A (t, θ), w A (t, θ), G M (t, θ), w M (t, θ), ω(t, θ), and λ(t, θ) of (t, θ) ∈ [0, ∞) × [−π, π], respectively.
In fact, this identification is possible in the following way. In general, an element x ∈ S corresponds one-to-one to an element θ ∈ [−π, π], i.e., so for any function f on S, there exists a corresponding periodic functionf on [−π, π] such that f (x) = f (x(θ)) =f (θ).
Then, integral of f on S is calculated by The distance between x = x(θ) and y = y(θ ) is calculated by where | · | abs stands for the absolute value.

Stationary solution
This section considers a homogeneous stationary solution and a linearized system of (18) on S with (19) 6 around the stationary solution.

Homogeneous stationary solution
We find the homogeneous stationary solution such that all variables Y , w A , G A , w M , G M , and ω are spatially uniform. Then, λ and φ should be .
From the first equation of (18), stationary Y becomes From the second equation of (18), stationary G A 6 In the following, we often refer to "(18) on S with (19)" simply as (18) for brevity when no confusion can arise.
Here, E α is defined by From the third equation of (18), stationary w A becomes Then, (22) together with (20) and (21) yield This gives, from (20), From the fourth equation of (18), stationary G M becomes Here E β is defined by From the sixth equation of (18) together with (21), (22) and (23), stationary ω is given by
Therefore, the sign of the eigenvalue (30) can be determined only by the sign of the content of the square brackets of (26).

Eigenvalue under extreme values of parameters
When so called "No black hole condition" holds, the eigenfunction can be stabilized for sufficiently high manufacturing transport cost (or sufficiently low preference for manufacturing variety). The no black hole condition same as Fujita et al. (2001, p.59) is given by the following theorem. See Appendix for proof.
Meanwhile, for sufficiently low manufacturing transport cost (or sufficiently high preference for manufacturing variety), the following theorem hold. See Appendix for proof.

Frequency dependence of eigenvalues
This section discusses how the eigenvalue (30) depends on the frequency of the eigenfunction (and the parameters of the model). In this section, the eigenvalues are computed numerically. In doing so, the range of the values of each parameter X(= τ A , τ M , η, or σ) is divided into an appropriate number N of equally spaced discrete points denoted by X 1 , X 2 , · · · , X N , where X i+1 − X i = ∆X for i = 1, 2, · · · , N − 1. Then, the eigenvalues are computed only over these discrete points. As a result, the illustrated curves in the following figures are those complemented by the drawing software. Figure 1 plots the eigenvalue (30) for several frequency numbers n with τ M on the horizontal axis and eigenvalues themselves on the vertical axis under µ = 0.5, η = 2.0, σ = 3.0, and τ A = 2.5. Here, ∆τ M ≈ 0.003 is set. Figure  1(a) shows that for n = 1 the eigenvalue is negative for any value of τ M . On the other hand, Figure 1(b) shows a more general situation for n = 2. That is, the eigenvalue is negative at sufficiently high values of τ M , positive at medium values of τ M , and negative again at low values of τ M . In this case, the eigenvalue (30) has two critical points τ * 1 M < τ * 2 M at which (30) changes its sign.

Critical curve
We investigate how the critical values depend on the parameters. For example, to examine the dependence on τ A , by plotting the points on the τ M -τ A plane where the value of the eigenvalue (30) is zero 8 , and we can draw a curve that divides the plane into stable and unstable areas. Here, the stable (unstable) area is the area where the homogeneous stationary solution is stable (unstable) for any pair of parameters (τ M , τ A ) in the region. Figure 2 shows the critical curve in τ M -τ A plane under µ = 0.5, η = 2.0, and σ = 3.0, and the frequency number is n = 3. Here, ∆τ A ≈ 0.0008 and ∆τ M ≈ 0.003 are set. It can be seen from Figure 2 that when τ A is sufficiently high, the two critical points τ * 1 M and τ * 2 M approach each other as τ A increases, and eventually merge and then disappear.  Figure 3 shows the critical curves plotted for different frequencies n under µ = 0.5, η = 2.0, and σ = 3.0. Here, ∆τ A ≈ 0.0008 and ∆τ M ≈ 0.003 are set. As can be seen from Figure 3a, the unstable area does not necessarily become larger as n increases 9 . However, if n is limited to even or odd numbers, we would observe that the unstable area expands or at least does not shrink as n increases (Figure 3(b) and 3(c)) 10 . 8 In the numerical computation, if the eigenvalue for the discrete point τM i is smaller (resp. greater) than that for τM i+1 , then we approximate τ * 1 M (resp. τ * 2 M ) by (τM i + τM i+1 )/2. 9 Clearly, the critical curves for n = 1 and n = 2 intersect. 10 See Remark below.  For the parameter η > 1, the critical curves can be drawn in the τ M -η plane. Figure 4 shows the critical curves on the τ M -η plane for some n under µ = 0.5, σ = 3.0, and τ A = 2.5. Here, ∆η ≈ 0.002 and ∆τ M ≈ 0.003 are set. Similar to the previous case in Figure 3, it can be seen from Figure 4 that when the values of η are sufficiently large, the two critical points τ * 1 M and τ * 2 M approach each other as η value becomes larger, and eventually merge and then disappear. The frequency dependence is also similar to the case of the critical curves on the τ M -τ A plane 11 .  For the parameter σ > 1, the critical curves under µ = 0.5, η = 2.0, and τ A = 2.5 are shown in Figure 5 12 . Here, ∆σ ≈ 0.03 and ∆τ M ≈ 0.01 are set. The shape of the critical curve looks quite different from previous cases in Figure 3 and 4. However, it is similar to the the previous cases qualitatively. Firstly, the two critical points move closer to each other as σ increases where the σ values are sufficiently large. Secondly, as the frequency number n increases only in the even or odd range, the two critical points move away from each other and the unstable region seems to expand or at least not to shrink 13 .  Remark: Numerical observations alone would not provide absolute confirmation of the claim that the unstable area is truly "monotonically expanding outward as the frequency number increases while remaining even or odd". This is not only because numerical observations are made only for a finite number of cases, but also because the critical curves corresponding to different frequencies can be so close to each other (especially on τ * 1 M side) that numerically they overlap exactly under some ranges of the parameters. For such cases, it would be impossible to infer the actual positional relationship of the critical curves, taking into account numerical errors.

Discussion and conclusion
We have considered a core-periphery model with differentiated agriculture and its transport cost on discrete and continuous multi-regional spaces. The stability of the homogeneous stationary solution of the model on a continuous periodic space has been investigated in detail. We have shown that the unstable eigenfunctions has two critical points at which the eigenvalue changes its sign. It has been observed numerically that the unstable area on a parameter plane between the two critical points generally tends to expand uniformly outward with the frequency increases only in even or odd numbers. Meanwhile, with the general increases in frequency, it has also observed that the unstable area does not necessarily expands uniformly outward.
Especially, what might be important in our results is the non-monotonic outward expansion of the unstable area with a mere increase in frequency n. In particular, the fact that the lower critical point τ * 1 M can move to the right with increasing n means that the economy can prefer to agglomerate into a relatively small number of regions rather than decentralize to a relatively large number of regions in a low range of the manufacturing transport cost. This is intuitively incompatible with the fact that the sufficiently low manufacturing transport cost makes the spatial structure more decentralized in the original two-regional model with differentiated agriculture (Fujita et al. (2001, p.110)). This non-trivial result can only be obtained by considering the continuous space (at least multi regional) model, and may reflect a richer spatial complexity than a simple two-regional relationship.
As noted in Remark to the previous section, our results rely on numerical observations. For a rigorous proof, it is necessary to compute analytically how the solution of the quadratic function as for H β n in (26) behaves with respect to changes in frequency n only within even numbers or only within odd numbers. However, there are difficulties that cannot be attributed to simple differential calculations.
When n is odd It is also immediate that b = 1 for α = 0 because α = 0 ⇔ H α n = 0, and Proof for Lemma 3 It is evident from (28) and Lemma 2.
Proof for Theorem 1 When the manufacturing transport cost goes to infinity, i.e., when β goes to infinity, the content in the square brackets of (26) becomes For there to be no black hole, (36) must be negative, for which it is sufficient that (1 − σ + µσ)(b − µ) < 0.
For α > 0, it is obvious that this limit is negative since B > 0 from Lemma 3.