Scientific Research

An Academic Publisher

Threshold Parameter for the Control of Unemployment in the Society: Mathematical Model and Analysis

**Author(s)**Leave a comment

1. Introduction

Securing a job or gain employment in our contemporary society is problematic for an average people and below; not to talk of desire job and underemployment is not something to rule out. In reference to Bureau of Labour Statistics defines unemployment as people who do not have a job, have actively looked for work in the past 4 weeks, and are currently available for work but could not get one or people who were temporarily laid off and were waiting to be called back to that job, Kimberly [1]. Galindro and Torress [2] referred to unemployment as a precarious social situation because a proportion of population normally struggles to maintain a minimum welfare and consumption level.

The unemployment is calculated as a percentage by dividing the number of unemployed individuals by all individual currently in labour force. During period of recession, an economy usually experiences a relatively high unemployment rate. According to the International Labour Organization (ILO), more than 200 million people globally or 6% of the work force were without job in 2012.

Economy of every nation is characterized by both active and inactive populations. The economically active ones are referred to as population willing and able to work and include those actively engaged in the production of goods and services and those who are employed, Muogbo and John-Akameli [3].

Nigeria has a challenge of accurate unemployment rates (i.e. data) to access. However according to Oyebade [4], Nigeria unemployment can be grouped into two categories; firstly, the older unemployed who lost their job through redundancy and bankruptcy and secondly, the younger unemployed, most of whom have never been employed. According to Awogbenle and Iwuamadi [5], the statistics from the man power board and the Federal Bureau of Statistics showed that Nigeria has a youth population of millions representing 60% of the total population of the country.

From Figure 1, it could be seen that, rate of unemployed youth has no year of declination in Nigeria despite efforts of Executives at various levels of government and in specific, 2017 N-POWER recruitment programme by President Muhammad Buhari administration. Also, 64 million of the youth are unemployed, while 1.6 million are under employed. The 1990-2000 data on youth unemployments showed that the largest group of unemployment is the secondary school leavers. Furthermore, 40% of the unemployment rate is among urban rural youth and 20 - 24 years old and 31% of the rate is among those aged 15 - 19 years old. Also the two third of the urban employed are rangeing from 15 - 24 years old. Also, the educated unemployed lie among young males with few dependents. There are relatively few secondary school graduates and the lower job

Figure 1. Trend of unemployment population (National Bureau of Statistics, Nigeria).

expectations of the primary school graduates. The causes of unemployment are heavily debated. Classical economics, new classical economics and the Austrian school of Economics argued that the market mechanisms are the reliable means of resolving means of unemployment; these theories against intervention imposed on the labour market from the outside such as unionization, bureaucratic work rules, minimum wages laws, taxes and other regulations that can claim discourage in the hiring of workers.

The problem of youth unemployment is very evident and severe in Nigeria. Every year, thousands of graduates are turn out for whom there are no jobs to accommodate them. Some are youth hawkers and bike riders who ordinarily would have found gained employment in some enterprises, or would have demonstrated their skills and resourcefulness if there are enabling environment and reliable management structure on ground. Instead, today youth have now shifted attention to cyber crime, terrorism and other criminal acts. It is the large number of youths who are unemployed capable of undermining democratic practices, constitute a serious threat if engaged by a political classes for clandestine activities.

In the study of youth unemployment activities especially in Nigeria, Adebayo [6], Alanana [7], Ayinde [8], Awogbenle and Iwuamadi [5], David and Ugocheckwu [9] have identified the main causes of youth unemployment in Nigeria. The rural-urban migration is one of the causes of unemployment as the rural area is no longer attractive, youth moves to urban areas with the motive of securing lucrative employment in the industries. In addition, is the misconception of social amenities in the urban centres. This mean that the rural areas are poorly consider in the allocation of social and economic opportunities.

Uddin and Osemengbe [10] also proposed a research that looked into the effects and solutions to youth unemployment in Nigeria. The data for their paper was collected from secondary sources using the descriptive approach of previous researchers and analysis of scholars to gather empirical data, their research revealed that unemployment in Nigeria among youths are caused by six major effects leading to communal clashes and the rise of such groups such as Boko Haram, Niger Delta Militant, Armed Robbery, Prostitution and child trafficking constituting hiccups to lifes and properties. Also, there research revealed that unemployment in Nigeria increases from 21.1% in 2010 to 23.9% in 2011 with youth unemployment at over 50% from 2011 to 2013, there is about increament of 16% unemployment growth rate in Nigeria. The impact of their paper also revealed that Government should invest on education heavily so that graduate will be self reliance instead of job seekers.

Another cause of unemployment in Nigeria is rapid growth rate. According tk worldometers, the current population of Nigeria is 196,364,552 as of Saturday 11th August, 2018, which was said to be the latest United Nations estimates while country meters read as of Sunday, 12th August, 2018, Nigeria population to be approximate 197.8 million. With this population, Nigeria is the most populous nation in Africa, the high population has resulted to rapid growth of labour force which is far outstripping the supply of jobs.

About mathematical contribution to knowledge of unemployment and method of analysis include; Pathan and Bhathawala [11] presented and analyzed a mathematical model for unemployment with fours dynamical variable; Unemployment $\left(u\left(t\right)\right)$ , Present job $\left(p\left(t\right)\right)$ , Vacancies $\left(v\left(t\right)\right)$ , and Employed person $\left(p\left(t\right)\right)$ . They analyzed the effect of the action of Government and private sector to control unemployment without and with the presence of delay, they also looked into the impact of self employment to decrease unemployment. Their analysis revealed the stability of non negative equilibrium point of the system, as well as numerical simulation for the validity of the analytical results.

Galindro and Torres [2] also proposes a simple mathematical model for unemployment. They worked with two compartmental models and the model was analyzed with real life data from Portugal. An optimal control was also formulated for the model and solved in which the results of their claim show realistic and non trivial solution.

Edogbanya and Ibrahim [12] Mathematical model of poverty and crime in human society. They developed a four compartmental models of poverty, crime, and incarceration using the system of non linear ordinary differential equation, they established the conditions that guarantee the stability of the model at crime free equilibrium and endemic equilibrium, they also determined the basic reproductive number R_{0} of the model and used it as a threshold condition, they also carried out the sensitivity analysis on the parameters and carry out numerical simulations with respect to the data.

Akinboro et al. [13] presented a paper on Numerical solution on SIR model using Differential transformation method and Variational iterative method. There studies investigated the application of Differential transformation method and variational iteration in finding the approximate solution of (SIR) model. SIR model are nonlinear system of Ordinary Differential Equation that has no analytic solution. The authors used general Lagrange multiplier concept to construct correction functional for VIM in solving the problem and tansformation of the original nonlinear system into numerical scheme called DTM. Their results revealed that both methods are in complete agreement, accurate and efficient for solving system of ODE. Many other authors have used VIM for model analysis and among are [14] [15] [16].

In view of all above, it is obvious that unemployment is one of the identified menace of contemprory society, especially, in Nigeria. To combat this problem, everybody is required to put positive effort and here I propose mathematical model to study the dynamics in population and set threshold to contain the challenge.

2. Mathematical Model Formulation and Analysis

2.1. Mathematical Model Formulation

The deterministic mathematical model was used to study the dynamics of unemployment. Individuals in the population are divided into compartments depending on the status of individual. The total human population at time t denoted by $P\left(t\right)$ , is divided into four sub-population of Unemployment $U\left(t\right)$ , and Unemployment compartment splitted into active, ${U}_{2}\left(t\right)$ and passive, ${U}_{1}\left(t\right)$ , $E\left(t\right)$ means Employment compartment, $R\left(t\right)$ denotes Retired compartment and lastly $V\left(t\right)$ represent vacancy/position available for recruitment compartment.

We assumed that there is an invisible interaction among all the human compartments and the vacancy compartment (this is represented with dotted line in the model diagram). The recruitment rate into the passive unemployment compartment is ${\pi}_{1}$ , ${\pi}_{2}$ is the recruitment rate into the active unemployment, ${k}_{1}$ is the rate at which individual moves from passive unemployment state, ${k}_{2}$ is the rate at which individual moves from active unemployment, to employment compartment, $\omega $ is the rate by which individual moves from employment class to retired. It is also assumed that people can still move from retired to both active and passive unemployment and the rate of their movements are denoted respectively by ${\beta}_{1}$ and ${\beta}_{2}$ . $\varphi $ is the vacancy creation of the employed people, ${\gamma}_{1}$ and ${\gamma}_{2}$ are the rate at which the employed people become unemployed probably, dissmisal resulting from misconduct or lack of productivity, $\delta $ are those that dies while working, $\eta $ is unpublicized available vacancies and $\mu $ is the natural death rate.

2.1.1. Flow Diagram of the Model

The model flow diagram (Figure 2) is described by the following system of nonlinear differential equations

$\frac{\text{d}{U}_{1}}{\text{d}t}={\pi}_{1}-\mu {U}_{1}-{k}_{1}{U}_{1}E+{\gamma}_{1}E+{\beta}_{1}R$ (1)

$\frac{\text{d}{U}_{2}}{\text{d}t}={\pi}_{2}-\mu {U}_{2}-{k}_{2}{U}_{2}E+{\gamma}_{2}E+{\beta}_{2}R$ (2)

$\frac{\text{d}E}{\text{d}t}=\left({k}_{1}{U}_{1}+{k}_{2}{U}_{2}\right)E-\left(\mu +\delta +{\gamma}_{1}+{\gamma}_{2}-\omega \right)E$ (3)

Figure 2. Unemployment Model Flow diagram.

$\frac{\text{d}R}{\text{d}t}=\omega E-\left(\mu +{\beta}_{1}+{\beta}_{2}\right)R$ (4)

$\frac{\text{d}V}{\text{d}t}=\varphi {U}_{2}+\gamma E-\eta V$ (5)

2.1.2. Assumptions of the Model

1) Only people of average and/or below (in a population) seek for job;

2) It is assumed that some unemployed individuals are not employable, yet, they seek for employment;

3) It is also part of assumption that some employed individuals have an establishment and also recruit from the pool of unemployment;

4) Some retired individuals added or return to unemployment population in either passive or active way;

5) For as long as people retire from job or loosing their job, there will always be vacancy and recruitment except it may not be publicized;

6) Available vacancy is somehow reserve for special people;

7) It is also assumed that the model considered able human being and dynamics of vacant positions; and

8) It is assumed that not everyone in the active compartment is looking for job.

2.2. Model Analysis

2.2.1. Basic Properties of the Model

The validity of any mathematical model depends whether the system of equations has a solution and communicate practical sense. If yes, is the solution unique? This subsection is concerned with examining if the system of equations has a solution(s) and possesses uniqueness, feasibility and positivity of solution properties.

Theorem 2.1.

Let $\Omega $ denote the region $0\le \alpha \le R$ , the system Equations (1)-(5) with initial conditions ${U}_{1}\left(0\right)>0$ , ${U}_{2}\left(0\right)>0$ , $E\left(0\right)>0$ , $R\left(0\right)>0$ and $V\left(0\right)\ge 0$ , exit, bounded, positively invariant and attracting for all $t>0$ .

Proof.

1) For existence and uniqueness of solution, we show that

$\frac{\partial {f}_{i}}{\partial {x}_{j}},i,j=1,2,3,4,5$ ,

are continuous and bounded in $\Omega $ .

From (1), ${F}_{1}={\pi}_{1}-\mu {U}_{1}-{k}_{1}{U}_{1}E+{\gamma}_{1}E+{\beta}_{1}R$

$\frac{\partial {f}_{1}}{\partial {U}_{1}}=-\mu -{k}_{1}E$ , $\left|\frac{\partial {f}_{1}}{\partial {U}_{1}}\right|=\left|-\mu -{k}_{1}E\right|<\infty $ , $\left|\frac{\partial {f}_{1}}{\partial {U}_{2}}\right|=\left|\frac{\partial {f}_{1}}{\partial V}\right|=0<\infty $ ,

$\left|\frac{\partial {f}_{1}}{\partial E}\right|=\left|-{k}_{1}{u}_{1}+{\gamma}_{1}\right|<\infty $ , $\left|\frac{\partial {f}_{1}}{\partial R}\right|=\left|{\beta}_{1}\right|<\infty $

From (2), ${F}_{2}={\pi}_{2}-\mu {U}_{2}-{k}_{2}{U}_{2}E+{\gamma}_{2}E+{\beta}_{2}R$

$\left|\frac{\partial {f}_{2}}{\partial {U}_{1}}\right|=0<\infty $ , $\left|\frac{\partial {f}_{2}}{\partial {U}_{2}}\right|=\left|-\mu -{k}_{2}E\right|<\infty $ , $\left|\frac{\partial {f}_{2}}{\partial E}\right|=\left|{\gamma}_{2}\right|<\infty $ ,

$\left|\frac{\partial {f}_{2}}{\partial R}\right|=\left|{\beta}_{2}\right|<\infty $ , $\left|\frac{\partial {f}_{2}}{\partial V}\right|=0<\infty $

Same approach for the remaining equations and consider

$\frac{\text{d}P\left(t\right)}{\text{d}t}={\pi}_{1}+{\pi}_{2}-\mu \left({U}_{1}+{U}_{2}+E+R\right)-\delta E$

$\frac{\text{d}P\left(t\right)}{\text{d}t}\le \pi -\mu P(t)$

where $\pi ={\pi}_{1}+{\pi}_{2}$

$\frac{\text{d}P}{\text{d}t}+\mu P\le \pi $ (6)

Integrating (6)

$P\left(t\right)\le {\text{e}}^{-\mu t}\left(\frac{\pi}{\mu}{\text{e}}^{\mu t}+C\right)$

Using initial condition (i.e. $P\left({t}_{0}\right)={P}_{0}$ at $t=0$ ), it implies that

$P\left(t\right)\le {P}_{0}{\text{e}}^{-\mu t}+\frac{\pi}{\mu}\left(1-{\text{e}}^{-\mu t}\right)$ (7)

$\Rightarrow {\mathrm{lim}}_{t\to \infty}P\left(t\right)=\frac{\pi}{\mu}$

Hence, $\frac{\pi}{\mu}$ is the upper bound of $P\left(t\right)$ .

Thus, when ${P}_{0}>\frac{\pi}{\mu}$ , the solution of $\left\{{U}_{1}\left(t\right),{U}_{2}\left(t\right),E\left(t\right),R\left(t\right)\right\}$ is bounded in the region of $\Omega $ .

Clearly, the partial derivatives of the whole system of equations exists, finite and bounded.

2) For feasibility region of solution

Let $P={U}_{1}+{U}_{2}+E+R\in {\mathbb{R}}_{+}^{4}$

$\frac{\text{d}P}{\text{d}t}=\frac{\text{d}{U}_{1}}{\text{d}t}+\frac{\text{d}{U}_{2}}{\text{d}t}+\frac{\text{d}E}{\text{d}t}+\frac{\text{d}R}{\text{d}t}\le \frac{{\pi}_{1}+{\pi}_{2}}{\mu}$

Therefore, $0\le P\le \frac{{\pi}_{1}+{\pi}_{2}}{\mu}$ and conclude that

$\Omega =\left\{\left({U}_{1}\left(t\right),{U}_{2}\left(t\right),E\left(t\right),R\left(t\right)\right)\subseteq {\mathbb{R}}_{+}^{4}:0\le P\le \frac{{\pi}_{1}+{\pi}_{2}}{\mu}\right\}$

is the feasible region.

3) For positivity of solution, from the proposition of the theorem, initial condition of the system is $\left\{\left({U}_{1}\left(0\right),{U}_{2}\left(0\right),E\left(0\right),R\left(0\right),V\left(0\right)\right)>0\right\}$ , then, we show that the solution set $\left({U}_{1}\left(t\right),{U}_{2}\left(t\right),E\left(t\right),R\left(t\right),V\left(t\right)\right)$ of the system remain non-negative for all $t>0$

$\frac{\text{d}{U}_{1}}{\text{d}t}={\pi}_{1}-{U}_{1}\left(\mu +{K}_{1}E\right)+{\gamma}_{1}E+{\beta}_{1}R$

$\frac{\text{d}{U}_{1}}{\text{d}t}\ge -\left(\mu +{k}_{1}E\right){U}_{1}$

$\int \frac{\text{d}{U}_{1}}{{U}_{1}}}\ge {\displaystyle \int -\left(\mu -{k}_{1}E\right)\text{d}t$

Integrating both side

${U}_{1}\left(t\right)\ge {U}_{1}\left(0\right){\text{e}}^{-\mu t}\ge 0$

Also,

$\frac{\text{d}{U}_{2}}{\text{d}t}={\pi}_{2}-\mu {U}_{2}-{k}_{2}{U}_{2}E+{\gamma}_{2}E+{\beta}_{2}R$

$\frac{\text{d}{U}_{2}}{\text{d}t}={\pi}_{2}-{U}_{2}\left(\mu +{k}_{2}E\right)+{\gamma}_{2}E+{\beta}_{2}R$

$\frac{\text{d}{U}_{2}}{\text{d}t}\ge -\left(\mu +{K}_{2}E\right){U}_{2}$

$\int \frac{\text{d}{U}_{2}}{{U}_{2}}}\ge {\displaystyle \int -\left(\mu +{k}_{2}E\right)\text{d}t$

$\mathrm{ln}{U}_{2}\ge -\left(\mu +{k}_{2}E\right)t+C$

${U}_{2}\ge {\text{e}}^{-\mu t}+C$

${U}_{2}\left(t\right)\ge {U}_{2}\left(0\right){\text{e}}^{-\mu t}\ge 0$

Similarly, $\frac{\text{d}R}{\text{d}t}=\omega E-\left(\mu +{\beta}_{1}+{\beta}_{2}\right)R$

$\frac{\text{d}R}{\text{d}t}\ge -\left(\mu +{\beta}_{1}+{\beta}_{2}\right)R$

$\int \frac{\text{d}R}{R}}\ge {\displaystyle \int -\left(\mu +{\beta}_{1}+{\beta}_{2}\right)R$

$R\left(t\right)\ge R\left(0\right){\text{e}}^{-\left(\mu +{\beta}_{1}+{\beta}_{2}\right)t}\ge 0$

By deducation,

$V\left(t\right)\ge V\left(0\right){\text{e}}^{-\eta t}>0$

Hence, for $t=0$ , all the state variables are non-negative. And this completes the proof.

2.2.2. Existence of Equilibria Points

We proceed to investigate the equilibrium points for the model, at equilibrium Equations (1)-(5) become

$0={\pi}_{1}-\mu {U}_{1}-{k}_{1}{U}_{1}E+{\gamma}_{1}E+{\beta}_{1}R$ (8)

$0={\pi}_{2}-\mu {U}_{2}-{k}_{2}{U}_{2}E+{\gamma}_{2}E+{\beta}_{2}R$ (9)

$0=\left({k}_{1}{U}_{1}+{k}_{2}{U}_{2}\right)E-\left(\mu +\delta +{\gamma}_{1}+{\gamma}_{2}+\omega \right)E$ (10)

$0=\omega E-\left(\mu +{\beta}_{1}+{\beta}_{2}\right)R$ (11)

$0=\varphi {U}_{2}+\gamma E-\eta V$ (12)

In the absent of recruitment, $E=0$ , $R=0$ , to have

${U}_{1}=\frac{{\pi}_{1}}{\mu},{U}_{2}=\frac{{\pi}_{2}}{\mu},V=\frac{\varphi {U}_{2}}{\eta}$

Therefore, the absent of recruitment equilibrium ( ${E}^{\text{*}}$ ) is obtained as

${E}^{\text{*}}=\left\{\left({U}_{1},{U}_{2},E,R,V\right):\left(\frac{{\pi}_{1}}{\mu},\frac{{\pi}_{2}}{\mu},0,0,\frac{\varphi {\pi}_{2}}{\eta \mu}\right)\right\}$ (13)

And persistent recruitment equilibrium, we assume $E=1$ , solving for the remaining variables, we get $R=\frac{\omega}{\mu +{\beta}_{1}+{\beta}_{2}}$ , ${U}_{1}=\frac{\left({\pi}_{1}+{\gamma}_{1}\right)\left(\mu +{\beta}_{1}+{\beta}_{2}\right)+{\beta}_{1}\omega}{\left({K}_{1}+\mu \right)\left(\mu +{\beta}_{1}+{\beta}_{2}\right)}$ , ${U}_{2}=\frac{\left({\pi}_{2}+{\gamma}_{2}\right)\left(\mu +{\beta}_{1}+{\beta}_{2}\right)+{\beta}_{1}\omega}{\left({K}_{2}+\mu \right)\left(\mu +{\beta}_{1}+{\beta}_{2}\right)}$ , 1, $\frac{\varphi {U}_{2}+\gamma}{\eta}$ .

Therefore the equilibrium of persistent recruitment denoted by ${E}^{\text{**}}$ obtained as

$\begin{array}{l}{E}^{\text{**}}=\{\left({U}_{1},{U}_{2},E,R,V\right):\frac{\left({\pi}_{1}+{\gamma}_{1}\right)\left(\mu +{\beta}_{1}+{\beta}_{2}\right)+{\beta}_{1}\omega}{\left({K}_{1}+\mu \right)\left(\mu +{\beta}_{1}+{\beta}_{2}\right)},\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\left({\pi}_{2}+{\gamma}_{2}\right)\left(\mu +{\beta}_{1}+{\beta}_{2}\right)+{\beta}_{1}\omega}{\left({K}_{2}+\mu \right)\left(\mu +{\beta}_{1}+{\beta}_{2}\right)},1,\frac{\omega}{\mu +{\beta}_{1}+{\beta}_{2}},\frac{\varphi {U}_{2}+\gamma}{\eta}\}\end{array}$ (14)

2.2.3. Threshold for the Recruitment

To examine the local stability of ${E}^{\text{*}}$ given by (13), is is imperative to obtain threshold parameter called basic reproduction number ( ${R}_{0}$ ) in the sense of epidemiology, which is defined as the average number of secondary infections caused by a typical infectious individual during its period of infectiousness in a completely susceptible population, (Samson et al., [17], Mastahun and Abdurahman, [18], Getachew et al., [19], Gabriel et al., [20] ). However, in this context, ( ${R}_{e}$ ) is referred to as the threshold for the recruitment into unemployed population.

We consider only the Employment and Retired compartments and with approach of next generation matrix (NGM)

$\frac{\text{d}}{\text{d}t}\left(\begin{array}{c}E\\ R\end{array}\right)=\left(\begin{array}{c}\left({k}_{1}{U}_{1}+{k}_{2}{U}_{2}\right)E\\ 0\end{array}\right)-\left(\begin{array}{c}\left(\mu +\delta +{\gamma}_{1}+{\gamma}_{2}+\omega \right)E\\ -\omega E+\left(\mu +{\beta}_{1}+{\beta}_{2}\right)R\end{array}\right)$

For which recruitment matrix (F) and transition matrix (V) respectively deduced as

$F=\left(\begin{array}{cc}{k}_{1}{U}_{1}+{k}_{2}{U}_{2}& 0\\ 0& 0\end{array}\right)$ and $V=\left(\begin{array}{cc}\mu +\delta +{\gamma}_{1}+{\gamma}_{2}+\omega & 0\\ -\omega & \mu +{\beta}_{1}+{\beta}_{2}\end{array}\right)$

to obtain threshold for recruitment, which is the spectral radius of the NGM, written mathematically as

${R}_{e}=\rho \left(F{V}^{-1}\right)=\frac{{k}_{1}{U}_{1}^{\ast}+{k}_{2}{U}_{2}^{\ast}}{\mu +\delta +{\gamma}_{1}+{\gamma}_{2}+\omega}$ (15)

${R}_{e}$ is the parameter that governs unemployment pool. When ${R}_{e}<1$ ; the situation for which the vacancy and employment population is more than unemployed individuals and it is closely observed. Otherwise, the economy of such soceity will greatly be affected.

2.2.4. Stability Analysis of the Model

Here, we carry out the steady states analysis; Absent of Recruitment steady state and the Persistent Recruitment steady state.

Therorem 2.2.

In the absence of recruitment equilibrium point, Equations (1)-(5) is locally stable if ${R}_{e}<1$ and unstable if otherwise.

Proof.

To check for the stability analysis of the recruitment free equilibrium of the model, we have to obtain the Jacobian matrix of the given Equations (1)-(5) as follow

$\mathbb{J}=\left(\begin{array}{ccccc}-\left(\mu +{k}_{1}E\right)& 0& -{k}_{1}{U}_{1}+{\gamma}_{1}& {\beta}_{1}& 0\\ 0& -\left(\mu +{k}_{2}E\right)& -{k}_{2}{U}_{2}+{\gamma}_{2}& {\beta}_{2}& 0\\ {k}_{1}& {k}_{2}& {k}_{1}{U}_{1}+{k}_{2}{U}_{2}-\left(\mu +\delta +{\gamma}_{1}+{\gamma}_{2}\right)& 0& 0\\ 0& 0& \omega & -\left(\mu +{\beta}_{1}+{\beta}_{2}\right)& 0\\ 0& \varphi & \gamma & 0& -\eta \end{array}\right)$ (16)

${\mathbb{J}}_{E=0}=\left(\begin{array}{ccccc}-\mu & 0& -{k}_{1}{U}_{1}+{\gamma}_{1}& {\beta}_{1}& 0\\ 0& -\mu & -{k}_{2}{U}_{2}+{\gamma}_{2}& {\beta}_{2}& 0\\ {k}_{1}& {k}_{2}& {k}_{1}{U}_{1}+{k}_{2}{U}_{2}-\left(\mu +\delta +{\gamma}_{1}+{\gamma}_{2}+\omega \right)& 0& 0\\ 0& 0& \omega & -\left(\mu +{\beta}_{1}+{\beta}_{2}\right)& 0\\ 0& \varphi & \gamma & 0& -\eta \end{array}\right)$ (17)

$\left|\mathbb{J}-\lambda \mathbb{I}\right|=\left|\begin{array}{ccccc}-\left(\mu +\lambda \right)& 0& -{k}_{1}{U}_{1}+{\gamma}_{1}& {\beta}_{1}& 0\\ 0& -\left(\mu +\lambda \right)& -{k}_{2}{U}_{2}+{\gamma}_{2}& {\beta}_{2}& 0\\ {k}_{1}& {k}_{2}& \left({k}_{1}{U}_{1}+{k}_{2}{U}_{2}-\left(\mu +\delta +{\gamma}_{1}+{\gamma}_{2}+\omega \right)\right)-\lambda & 0& 0\\ 0& 0& \omega & -\left(\mu +{\beta}_{1}+{\beta}_{2}+\lambda \right)& 0\\ 0& \varphi & \gamma & 0& -\left(\eta +\lambda \right)\end{array}\right|$

$\begin{array}{c}\left|J-I\lambda \right|={\left(\mu +\lambda \right)}^{2}\left[\left({k}_{1}{U}_{1}+{k}_{2}{U}_{2}-\left(\mu +\delta +{\gamma}_{1}+{\gamma}_{2}+\omega \right)\right)-\lambda \right]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\times \left(\mu +{\beta}_{1}+{\beta}_{2}+\lambda \right)\left(\eta +\lambda \right)\\ =0\end{array}$

Solving the above expression to obtain values for the λ’s, we have;

$\begin{array}{l}{\lambda}_{1}=-\mu ,\text{\hspace{0.17em}}{\lambda}_{2}=-\mu ,\text{\hspace{0.17em}}{\lambda}_{3}={k}_{1}{U}_{1}^{\ast}+{k}_{2}{U}_{2}^{\ast}-\left(\mu +\delta +{\gamma}_{1}+{\gamma}_{2}+\omega \right),\\ {\lambda}_{4}=-\left(\mu +{\beta}_{1}+{\beta}_{2}\right),\text{\hspace{0.17em}}{\lambda}_{5}=-\eta \end{array}$ (18)

Condition for stability is that, all eigenvalues of the Jacobian matrix must be less than zero [20]. Hence, ${\lambda}_{3}$ require to be less than zero (i.e. removal rate to employment must be greater than recruitment rate to unemployment) for stability.

Conversely, if ${\lambda}_{3}>0$ , then, the case employment free equilibrium will be unstable and further analysis would be required, which it not the interest of this research.

2.2.5. Application of Variational Iteration Method (VIM)

Applying the Semi-analytic Numerical Scheme; VIM and employ concept used in [13] for the correctional functionals on our model as:

$\begin{array}{l}{U}_{1,K+1}\left(t\right)={U}_{1}\left(K\right)+{\displaystyle {\int}_{0}^{t}{\lambda}_{1}\left(t\right)}(\frac{\text{d}{U}_{1}\left(K\right)}{\text{d}t}-{\pi}_{1}+\left(\mu +{k}_{1}E\left(K\right)\right){U}_{1}\left(K\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\begin{array}{c}\\ \end{array}-{\gamma}_{1}E\left(K\right)+{\beta}_{1}R\left(K\right))\text{d}t\\ {U}_{2,K+1}\left(t\right)={U}_{2}\left(K\right)+{\displaystyle {\int}_{0}^{t}{\lambda}_{2}\left(t\right)}(\frac{\text{d}{U}_{2}\left(K\right)}{\text{d}t}-{\pi}_{2}+\left(\mu +{k}_{2}E\left(K\right)\right){U}_{2}\left(K\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\begin{array}{c}\\ \end{array}-{\gamma}_{2}E\left(K\right)+{\beta}_{2}R\left(K\right))\text{d}t\end{array}$

$\begin{array}{l}{E}_{K+1}\left(t\right)=E\left(K\right)+{\displaystyle {\int}_{0}^{t}{\lambda}_{3}\left(t\right)}(\frac{\text{d}E\left(K\right)}{\text{d}t}-\left({k}_{1}{U}_{1}\left(K\right)+{k}_{2}{U}_{2}\left(K\right)\right)E\left(K\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\begin{array}{c}\\ \end{array}+\left(\mu +\delta +{\gamma}_{1}+{\gamma}_{2}\right)E\left(K\right)+{\beta}_{2}R\left(K\right))\text{d}t\\ {R}_{k+1}\left(t\right)=R\left(K\right)+{\displaystyle {\int}_{0}^{t}{\lambda}_{4}\left(t\right)}\left(\frac{\text{d}R\left(K\right)}{\text{d}t}-wE\left(K\right)+\left(\mu +{\beta}_{1}+{\beta}_{2}\right)R\left(K\right)\right)\text{d}t\\ {V}_{K+1}\left(t\right)=V\left(K\right)+{\displaystyle {\int}_{0}^{t}{\lambda}_{5}\left(t\right)}\left(\frac{\text{d}V\left(K\right)}{\text{d}t}-\varphi {U}_{2}\left(K\right)-\gamma E\left(K\right)+\eta V\left(K\right)\right)\text{d}t\end{array}$ (19)

it is shown from Akinboro et al., [13] that ${\lambda}_{1}\left(t\right)=-1$ , ${\lambda}_{2}\left(t\right)=-1$ , ${\lambda}_{3}\left(t\right)=-1$ , ${\lambda}_{4}\left(t\right)=-1$ , ${\lambda}_{4}\left(t\right)=-1$ , substituting the multipliers into the correctional functional above, then, we have the following iterative formula

$\begin{array}{l}{U}_{1,K+1}\left(t\right)={U}_{1}\left(K\right)-{\displaystyle {\int}_{0}^{t}(\frac{\text{d}{U}_{1}\left(K\right)}{\text{d}t}-{\pi}_{1}+\left(\mu +{k}_{1}E\left(K\right)\right){U}_{1}\left(K\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\begin{array}{c}\\ \end{array}-{\gamma}_{1}E\left(K\right)+{\beta}_{1}R\left(K\right))\text{d}t\end{array}$

$\begin{array}{l}{U}_{2,K+1}\left(t\right)={U}_{2}\left(K\right)-{\displaystyle {\int}_{0}^{t}(\frac{\text{d}{U}_{2}\left(K\right)}{\text{d}t}-{\pi}_{2}+\left(\mu +{k}_{2}E\left(K\right)\right){U}_{2}\left(K\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\begin{array}{c}\\ \end{array}-{\gamma}_{2}E\left(K\right)+{\beta}_{2}R\left(K\right))\text{d}t\end{array}$

$\begin{array}{l}{E}_{K+1}\left(t\right)=E\left(K\right)-{\displaystyle {\int}_{0}^{t}(\frac{\text{d}E\left(K\right)}{\text{d}t}-\left({k}_{1}{U}_{1}\left(K\right)+{k}_{2}{U}_{2}\left(K\right)\right)E\left(K\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\begin{array}{c}\\ \end{array}+\left(\mu +\delta +{\gamma}_{1}+{\gamma}_{2}\right)E\left(K\right)+{\beta}_{2}R\left(K\right))\text{d}t\\ {R}_{k+1}\left(t\right)=R\left(K\right)-{\displaystyle {\int}_{0}^{t}\left(\frac{\text{d}R\left(K\right)}{\text{d}t}-wE\left(K\right)+\left(\mu +{\beta}_{1}+{\beta}_{2}\right)R\left(K\right)\right)\text{d}t}\\ {V}_{K+1}\left(t\right)=V\left(K\right)-{\displaystyle {\int}_{0}^{t}\left(\frac{\text{d}V\left(K\right)}{\text{d}t}-\varphi {U}_{2}\left(K\right)-\gamma E\left(K\right)+\eta V\left(K\right)\right)\text{d}t}\end{array}$ (20)

3. Numerical Simulation and Discussion of Results

3.1. Numerical Simulation

Implementing (20) on maple 18 with initial conditions ${U}_{1}\left(0\right),{U}_{2}\left(0\right),E\left(0\right),R\left(0\right)$ and $V\left(0\right)$ and values in Table 1, the following infinite series solutions were obtained for K = 4.

$\begin{array}{l}{U}_{1}\left(t\right)=60-1186.0t-20248.83{t}^{2}+196082.19{t}^{3}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+12000399.09{t}^{4}+74310736.32{t}^{5}+\cdots \\ {U}_{2}\left(t\right)=40-524.50t-10934.45{t}^{2}+9122.34{t}^{3}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+4435444.89{t}^{4}+58674358.19{t}^{5}+\cdots \end{array}$

$\begin{array}{l}E\left(t\right)=35+1838.60t+32191.71{t}^{2}-193491.25{t}^{3}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}-16490718.82{t}^{4}-136637622.6{t}^{5}+\cdots \\ R\left(t\right)=5+27.50t+901.42{t}^{2}+7753.30{t}^{3}-60749.11{t}^{4}-1798174.31{t}^{5}+\cdots \\ V\left(t\right)=10+56.00t+568.01{t}^{2}+4747.13{t}^{3}-42286.075{t}^{4}-1147347.94{t}^{5}+\cdots \end{array}$ (21)

3.2. Discussion of Results

This article recorded the following findings:

・ When there is no employment, it is practically sensible that there would be no retirement but vacancy exists due to active unemployment and may/maynot be publicize, Equation (13).

・ As employment exist, the result obtained illustrates retirement, vacancies and population of unemployment (both passive and active), all exist, Equation (14).

・ A threshold that is to establsih stability of unemployment was obtained and required to be well manage, Equation (15). When ${R}_{e}<1$ , implies people getting job are more than those without job, otherwise, the economy of such community will face hardish.

・ It could be seen that all eigenvalues of the characteristics equation give expression less than zero except ${\lambda}_{3}$ , which for stability, must be less than zero and this said ${\lambda}_{3}$ is the threshold (implies the assurance of resolving problem of unemployment in contemporal world depends on the threshold) obtained for unemployment, see Equation (18). This is in support of result obtained is (15) and it importance.

・ VIM was employed to obtain an approximate solution of the propose model and the solution was hypothetically fitted with values to produce results in Table 2, which shows how the population behaves as time goes by. Starting from year 1 to year 5. At first year, it showed that the available vacancies is not in proportion with the available population and it is practically sensible that, out of the total population of unemployment in the first year, almost half of the population are get employed, which showed that the remaining unemployed individuals added up to the next year population and that is why the population of unemployment keep rising as time goes.

4. Summary, Conclusion and Recommendations

4.1. Summary

Here is summary of the research carried out on mathematical modeling of unemployment with retirement and conclusion made based on our findings.

Table 1. Definition of parameters and Simulation value.

Table 2. The model solution obtained by V.I.M. was arbitrarily truncated and hypothetically fitted to observe the population of the dynamical model.

In this research, we are able to modify the model of [2] features the dynamics of unemployment with consideration of retirement. We established the existence and uniqueness of solution of the model. The assurance of resolving problem of unemployment in the population was also, achieved. We further obtained the threshold that must be met in order for stability of the unemployment population. The last is the use of VIM to validate the analytic solutions of the model. A radical behaviour of the model was observed if quick intervention is not surface.

4.2. Conclusions

The cooperate environment is now very dependent on skills and technical know-how, western education and certificate(s) is/are no longer enough/unenthusiastic to provide food on the table. Hence, this article concluded that:

1) Every individual have contribution to end unemployment by developing self skills of entrepreneur or learn one (this can be view from Equation (13)).

2) Important threshold (see Equation (15)) to curtail challenges of unemployment was established and beseeched government and policy makers to introduce policy/strategy in proportion to exponential population growth.

3) It is imperative for the employed individuals to prepare for retirement as pension scheme is not sufficient to cater for life after 35years of service (Equation (18)).

4.3. Suggestion for Further Research

This article studied dynamics in mathematical sense, unemployment and established important parameter, called threshold, unemployment is said to be stable for ${R}_{e}<1$ (when no recruitment except active unemployed created job) and ${R}_{e}>1$ will also stable when there is employment (i.e. when everybody put effort and create job) but we did not consider the case ${R}_{e}=1$ , which could be refered to as bifurcation analysis in the mathematical view (when total number of vacant position is equal to number of population). Also, the model can be modified based on asumptions different from those here considered.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

*Journal of Applied Mathematics and Physics*,

**6**, 2563-2578. doi: 10.4236/jamp.2018.612214.

[1] |
Kimberly, A. (2018) Unemployment, Its Causes and Consequences. https://www.thebalance.com/what-is-unemployment-3306222 |

[2] | Galindro, A. and Torress, D.F.M (2017) A Simple Mathematical Model for Unemployment: A Case Study & Portugal with Optimal Control. Statistics, Optimization and Information Computing, 6, 116-129. |

[3] | Muogbo, U.S. and John-Akamelu, C.R. (2018) Impact of Enterpremeural Skills in Reducing Youth Unemployment in Nigeria. European Journal of Business, Economics and Accountancy, 6, 1-12. |

[4] | Oyebade, S.A. (2003) Education and Unemployment of Youth in Nigeria: Causes, Impact and Suggestions. National Economics Empowerment Development Strategy (NEEDS), 45-65. |

[5] | Awogbenla, A.C. and Iwuamadi, K.C. (2010) Youth Unemployment: Enterprenuership Development Programe as an Intervention Mechanisms. African Journal of Business Management, 4, 831-835. |

[6] | Adebayo, A. (1999) Youth Unemployment and National Directorate of Employment, Self Employmet Programme. National Journal of Economics and Social Studies, 41, 81-102. |

[7] |
Alanana, O.O. (2003) Youth Unemployment in Nigeria: Some Implications for the Third Millenium. Global Journal of Social Sciences, 2, 21-26. https://doi.org/10.4314/gjss.v2i1.22763 |

[8] | Ayinde, O.E. (2008) Empirical Analysis of Agricultural Growth and Unemployment in Nigeria. Africa Journal of Agricultural Research, 3, 465-468. |

[9] | David, I. And Ugockukwu, M.U. (2018) Youth Employment in Nigeria and Impregnable but Artificial Walls: The Urgency of a New and Inclusive Country. International Sociological Association. XIX World Congress of Sociology, Toronto, 15-21 July 2018, 1-21. |

[10] | Uddin, P.S.O. and Osemengbe, O. (2013) The Causes, Effect and Solution to Youth Unemployment Problems in Nigeria. Journals of Emerging Trends in Economics and Management Sciences (JETEMS), 4, 397-402. |

[11] | Pathan, G. and Bhathawala, P.H. (2017) Mathematical Model of Unemployment—An Analysis with Delay. Global Journal of Mathematical Sciences: Theory and Title in Practical, 9, 225-237. |

[12] | Edogbanya, H. and Ibrahim, M.O. (2016) Mathematical Model of Poverty and Crime in Human Society. Unpublished Ph.D. Thesis, Department of Mathematics, Faculty of Physical Sciences, University of Ilorin, Ilorin. |

[13] | Akinboro, F.S., Alao, S. and Akinpelu, F.O. (2014) Numerical Solution of SIR Model Using Differential Transformation Method and Variational Iteration Method. General Mathematics Notes, 22, 82-92. |

[14] |
Rafei, M., Daniali, H. and Ganji, D.D. (2007) Variational Iteration Method for Solving the Epidemic Model and the Prey Predator Problem. Journals on Applied Mathematics and Computation, 186, 1701-1709. https://www.sciencedirect.com/ https://doi.org/10.1016/j.amc.2006.08.077 |

[15] | He, J.H. (2000) Variational Iterative Method for Autonomous Ordinary Differential Systems. Application of Mathematics and Computation, 114, 115-123. |

[16] |
Abbasbandy, S. and Shivanian, E. (2009) Application of the Variational Iteration Method for System of Non Linear Volterra Integro Differential Equation. Mathematical and Computer Application, 14, 147-158. https://doi.org/10.3390/mca14020147 |

[17] | Samson, O., Maruf, A.L. and Olawale, S.O. (2016) Stability and Sensitivity Analysis of a Deterministic Epidemiological Model with Pseudo-Recovery. IAENG International Journal of Applied Mathematics, 46, 2. |

[18] |
Mamatjan, M. and Xamxinur, A. (2017) Optimal Control of an HIV/ADIS Epidemic Model with Infective Immigration and Behavioural Change. Applied Mathematics, 8, 87-105. https://doi.org/10.4236/am.2017.81008 |

[19] | Getachew, T.T., Oluwole, D.M. and David, M. (2017) Modelling and Optimal Control of Typhoid Fever Disease with Cost-Effectiveness Strategies. Computational and Mathematical Methods in Medicine, Hindawi, 1-16. |

[20] | Gabriel, O., Joseph, K.K. and John, M.M. (2016) Cost Effectiveness Analysis of Optimal Malaria Control Strategies in Kenya. Moi University, Kenya. |

[21] |
Nigeria Youth Unemployment. Nigeria Bureau of Statistics. https://tradingeconomics.com/nigeria/youth-unemployment-rate |

Copyright © 2019 by authors and Scientific Research Publishing Inc.

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.