A simple mathematical model for unemployment: a case study in Portugal with optimal control

We propose a simple mathematical model for unemployment. Despite its simpleness, we claim that the model is more realistic and useful than recent models available in the literature. A case study with real data from Portugal supports our claim. An optimal control problem is formulated and solved, which provides some non-trivial and interesting conclusions.


Introduction
Unemployment is an extremely severe social and economic problem, born from the difference between demand and supply of the labour market and sometimes emphasized by population's growth. The unemployed population can be defined as the portion of able citizens that desire to work (known as the active population) but, unfortunately, due to insufficient supply, are deprived from working. Generally, unemployment is a precarious social situation since a portion of the population normally struggles to maintain a minimum welfare and consumption level. Simultaneously, regarding the macroeconomic perspective, higher unemployment rates intensify the pressure on social protection measures, e.g., unemployment subsidies, and, consequently, the associated government expenditure. There are numerous policies that interact directly with a country's level of unemployment or, symmetrically, with the level of employment offered to the citizens. As an example, we have the establishment of a minimum wage and currency devaluation. Here, we propose and investigate a simple new model that describes well the current labour market in Portugal.
Regarding available literature focused on the subject, we go back to 2011, time when Misra and Singh [5] proposed a non-linear mathematical model of unemployment based on Nikolopoulos and Tzanetis previous work of 2003 [8]. More precisely, they suggested a model for housing allocation of a homeless population due to a natural disaster, described by a system of ordinary differential equations with the following three variables: unemployed population, employed individuals, and temporarily workers [5]. They analyse the equilibrium of the model, using the stability theory for differential equations, and perform a few numerical simulations, concluding that the unemployment battle may need immediate measures: they predict that the unemployment rate may rise quickly and, if those high unemployment values are reached, then it might be very difficult to overcome that much bigger problem in the future [5]. Misra and Singh developed further their empirical work, and in 2013 they

Data analysis and the Munoli and Gani 2016 model
As part of our goal, we collected data from IEFP (Instituto do Emprego e Formação Profissional), a Portuguese government institution focused on providing education and support for the unemployed population [12] and Bank of Portugal [11]. We compiled the number of unemployed persons (U ), the unemployment rate (U R) and vacancies available at the IEFP (D), concerning the time period from January 2004 until June 2016 with monthly frequency. A total of 150 observations from each variable were collected, U t , U R t , D t , t = 1, . . . , 150, and three new derived variables were defined: • unemployment change rate • employment change rate Formulas (2) and (3) indicate the rate of change of unemployed and employed people, respectively. We could not obtain the total number of employees in Portugal during our time frame. Therefore, we proceeded with an indirect calculation (1) using the variables that we actually could gather: the unemployment rate U R and the total number of unemployed people U .
3 Munoli and Gani (2016) [7] present a model that tries to emulate an unemployment environment. Their model consists of three ordinary differential equations, considering the unemployed (U ), the employed (E), and available vacancies (V ): The parameters and variables of model (4) are described on Table 1. At Munoli and Gani (2016) [7], the initial Denotes the diminution rate of vacancies due to lack of funds conditions were given by U (0) = 10000, E(0) = 1000 and V (0) = 100. A time line of 150 units (t = 150) was considered [7]. Values of the other parameters were fulfilled according with Table 2. We replaced the initial  [7] conditions suggested by the authors of [7] by the ones given by the real data from Portugal. Precisely, we fixed the initial values to the ones of Portugal at January 2004, when the number of unemployed people was 464450 (U ), employed was 6450694 (E, according with (1)) and the available vacancies at the time were 4848 (V ). We can easily state that the suggested model from Munoli and Gani (2016) [7] does not replicate properly the real data from Portugal. Regarding the employment and unemployment (see Figures 1a and 1b, respectively) the model of [7] kind of implode, since the unemployed/employment values dramatically drop until the end of the time period. Munoli and Gani (2016) [7] differential system (4) also suggests an exceptional increase on the supply of job vacancies with a smoother decrease until the end of the time period, a statement that is not supported by the Portuguese real data. Actually, the real number of vacancies presents a smooth fluctuation and a shy tendency to increase over time (see Figure 1c).

New unemployment model and simulations
Given the weak results of Section 2, our main goal is to create a new mathematical model that explains more accurately the unemployment environment. With this in mind, we proceed with a few changes in Munoli and Gani (2016) [7] model to achieve the desired result. First of all, we consider that the number of vacancies should be an exogenous variable and not given by any specific differential equation, as stated on previous works. The inherent fluctuation and apparent pattern lead us to fit our data in order to achieve a trustworthy representation of this variable. We fit our data with a 3rd degree Fourier function (see Appendix A), obtaining a reasonable goodness of fit (precisely, a R-square of 0.8046, as given in Appendix A). Figure 2 shows the fitted function and the actual vacancies data.
Having in mind the results of Sîrghi et al. of 2014 [10], we decided to include one variable that compiles the employment created due to the increase of the unemployment. Bigger unemployment rates signals employers to hire at a lower wage and, as a consequence, they are also able to hire more workers. We computed the p-values for Pearson's correlation using a Student's t-distribution for a transformation of the correlation using the corr(X,Y) command of MatLab. This function is exact when X and Y are normal. To retrieve such correlation value, we used the two transformed variables RCU t and RCE t (formulas (2) and (3)), obtaining the value of 0.7161. We also found that the constant rate, at which the number of unemployed persons is increasing continuously, as assumed by Munoli and Gani (2016) [7] and Misra and Singh (2013) [6], is quite small regarding the Portuguese population. For this reason, we have increased the value to 90000. We also note that the other crucial point explaining the under-achievement of previous models with respect to Portuguese data (implosion and general shrinkness of the population, see Section 2), is the absence of a value of a constant rate at which the number of employed persons is increasing continuously in order to recoup the loss of people within the system. Considering the stated above, we propose here the following model for unemployment, described by a system of two ordinary differential equations: The meaning of the variables and parameters of our model (5) is given in Table 3. The values used in our simulations are given in Table 4. Denotes the diminution rate of vacancies due to lack of funds ρ Rate of employment increase due to labour-force wage devaluation The feasible region of model (5) is given by Next result gives the positivity invariance of the feasible region (6).

Theorem 1
The set Ω defined by (6) is a region of attraction for the model system (5).
Given model (5) and the parameter values of Table 4, we carried out the simulation in Matlab (see Appendix A). Since the required T to smooth our differential equation was 81 observations, we needed to compress our observed real values (150 observations) in order to achieve a graphical comparison (see Figure 3). As we can see, our model fits the observed data much better than the results obtained from previous models available in the literature. Even though there are a few opposite fluctuations, our model achieves a much more steady environment.

Equilibrium analysis
To compute the equilibrium points of the proposed model (5), one needs to solve the following system: Direct calculations show that there exists one equilibrium point E b = (U * , E * ) only, given by Remark 1 All the parameters that appear in (7) are strictly positive. Therefore, the numerators of U * and E * are always positive. The only possibility for U * and E * to be negative would be to have which is not a reasonable scenario since the variable V , representing the available vacancies, is way bigger than all other parameters appearing in the denominators of U * and E * , reasonably valued in the interval [0, 1].

Stability analysis
We now study the local stability of the equilibrium E b found in Section 3.1. To achieve this, we compute the so called variational matrix M of our designated model (5): The characteristic equation associated with our 2 × 2 matrix (8) is The Routh-Hurwitz criterion for second degree polynomials asserts that (9) has all the roots in the left half plane (and the system is stable) if and only if both coefficients (10) are positive [2]. We just proved the following result.

Remark 2
For the values used to describe the Portuguese reality of unemployment, one has a 2 = 0.00000090 V + 0.0033239, which is strictly positive because V is always non-negative. It follows that the unique equilibrium point of our system is locally asymptotically stable.

Optimal control
Portugal is a country with serious unemployment issues during the last decade and the government of Portugal was forced to apply intervention policies in this particular area. Regarding this subject, the adoption of internships and hiring support measures (policies where the government contributes with a share of the worker's salary during a pre-established period, normally one year) have been in force since 1997, with variable magnitude until nowadays. Facing more severe unemployment problems at the 2007/2008 crisis, those measures became quite popular as a fundamental axle regarding the battle against unemployment and integration in the labour market of the recent graduates. With reference to the bibliography on this subject, there are a few empirical works that try to answer or explain the following difficult question: "Does the supply of internships fight the long-term unemployment?" The work of Silva et al. (2016) [9] focus on the impact of the internships in the unemployment of graduates compared to other age-similar groups. Another study, Barnwell (2016) [1], addresses the effectiveness of the internship component in the increasing employability of graduates. However, we are not aware of any empirical work that responds concretely to the aforementioned question. Using our representative model of the labor market reality, we now introduce two controls into (5): The first control function, u 1 , refers to the unitary supply of internships or support measures; while the second control function, u 2 , represents other alternative indirect measures such as lowering corporate tax rates. The offer of an internship, represented by the control measure u 1 , has a direct or immediate impact based on the simple premise that an unemployed worker shifts to the employed group due to this incentive. This variable is scaled between −40000 and 40000 because it is possible to add or withdraw internships already operating in the market. The cost of each internship is registered by the monetary value of the support plus all the administrative costs inherent to the planning and execution of the internship. The magnitude of indirect benefits, denoted by the control variable u 2 , interacts directly with the exogenous function since those measures affect the natural creation of employment. We settled its value between 0 and 1, being its cost interpreted in the monetary unit (m.u.) of internships. The optimal control problem is as follows: 12) subject to the constraints on the control values the initial conditions the terminal conditions 5000000 ≤ U (150) + E(150) ≤ 8000000 (15) and the state constraint Since we know that the supply/withdraw of a unitary control u 1 represents a new employed/unemployed in the system, and a government financial cost/gain inherent to this measure, we set B as the reference value equal to 1.

9
In order to settle the value of A, we establish an ideal ratio of 20 to 1, stating that an unemployed person has the similar cost as 19 new employees, using 5% as a target of the utopian level of unemployment. The value of C is set at 40000, representative of 40000 m.u. expressed in internship currency so that when u 2 is equal to 1 (maximum value) we are stating that we are investing the value/cost of 40000 internships in indirect measures. The constraint (16) keeps the unemployment rate below 12% while (15) assures a reasonable level of active labour force between 5 and 8 million people. We have chosen the value 12% since this is the minimum unemployment rate in the last 5 years of our real data. The initial conditions (14) of employment/unemployment level agree with collected data. The optimal control problem (12)-(16) is far from trivial and we were not able to find its analytical solution through application of the Pontryagin Maximum Principle for optimal control problems with state constraints. Instead, we solved the optimal control problem (12)-(16) numerically using the ACADO Toolkit [4] -see Appendix B. The results are given in Figures 4 and 5. As we see, looking to the graphics in Figure 4, the model suggests a moderate (about 0.5) adoption of the u 2 control from the beginning of the period until t = 50, point at which the number of unemployed people approaches its minimum in our study. At this point it is suggested to switch the selected control: the u 2 control shrinks to zero while a substantial enlarged policy of internships (u 1 ) is suggested during the time-frame from t = 50 until approximately t = 110. The reason for this policy, during the period t ∈ [70, 110], might be related to the employment minimum value. Finally, from t = 110 until the end of the simulation, we assist to a new rise in unemployment levels. The optimal control approach points out a brief massive (u 2 = 1) followed by a moderate new supply of indirect incentives u 2 (approximately 0.2) and an internship total contraction, that is, u 1 = −40000. The total cost of applying the controls during our time-frame suggests a slight increase in the total investment up to t = 110 culminating in a final saving and financial recovery. We note that the application of controls slightly tip the expenditure in approximately 1000 more internships (monthly) during the 150 month study period. Finally, comparing the actual data with our optimal control simulation, the level of unemployment regarding the simulation surpasses slightly the real data from the beginning of the simulation until t = 60. During the remaining period, the optimal control approach avoids the severe unemployment period from the actual Portuguese data between t = 100 and t = 110, which scores values above 17%. Overall, and considering the mean over the analysed period, the simulation with optimal controls scores 9.8% slightly better than the real data value of 11.5%, as we see in Figure 5.

Conclusion
In this work, we propose a simple mathematical model that describes more accurately the real data of unemployment from Portugal in the period 2004-2016. At first, we considered the proposed model from Munoli and Gani (2016) [7], fitting their initial conditions with our data. We found out that the model suggested by Munoli and Gani in 2016 does not describes at all the Portuguese unemployment environment. Because of this, we performed a few changes to the Munoli and Gani model [7], eliminating one of the differential equations, adding an unemployment/wage correlation feature, and adjusting the natural rate of unemployed/employed people in order to control and balance the inner population of the model. Simulations of our model show a much more accurate emulation of the Portuguese unemployment environment. Our results also show a slight decrease on the overall labour force (unemployed plus employed) and, since the Portuguese total population smoothly increased over the last decade, that might signals an higher percentage of inactive population, which may underpin pressure on more social protection measures. Other reading may consist on the premise that the number of unknown unemployed people is increasing and standing apart from the government official records. From the application of optimal control, we can state the following interesting conclusions: the indirect policies should be the predominant method of avoiding unemployment, whereas the supply of internships should be the main choice when the total level of employment offered is low. A possible explanation why to avoid internships in high unemployment periods might be correlated to the severe labour devaluation (considered when we developed the considered model) during these seasons, since the available labour force is already cheap enough and the internships supply might be offering jobs with an increased government expenditure, due to the whole administrative costs, to people that might turned out    to be employed anyway. Moreover, our study also suggests that governmental policies should be performed mainly during favourable periods in order to avoid high unemployment levels in future crisis. As future work, we plan to apply different methods of optimal control to our model in order to seek even more solid ways and tools to control the unemployment issue. For that, we may consider to split the unemployment class into two: unemployed that currently receive welfare from government; and unemployed which do not receive any financial support. These two classes present different government concerns: the first one emphasizes financial pressure and the second one social and well-being pressure. Another interesting study might be to include nonactive population, like retired people, and study the optimal control regarding the social security financial health.  We obtained a 1737 × 3 matrix (denoted above by xa). Thus, we formatted the space in order to make the comparison with our data (a 150 × 3 matrix): Using the MatLab fitting tool, we find a 3rd degree Fourier function that fits the total vacancies real data quite well: General model Fourier3: f(x) = a0 + a1 * cos(x * w) + b1 * sin(x * w) + a2 * cos(2 * x * w) + b2 * sin(2 * x * w) + a3 * cos(3 * x * w) + b3 * sin (