The health workforce is an important factor in the functioning of the health system. Therefore, to have an efficient health system, sufficient workforce, and optimal ratio, it is necessary to ensure that practitioners are motivated enough to perform their assigned functions well (1). Human capital is also a key component of sustainable development (2), and physicians as human resources and human capital are among the most important factors affecting health system performance and sustainable growth and development. Accordingly, the imbalance between the supply and demand of the physician leads to the inefficiency of the health system. Shortage of physicians leads to major problems in accessing health services, and decreasing health status and physician surplus lead to the high expenditures of public health and waste of resources (3). Therefore, careful planning is necessary for the human capital of health, particularly providing an adequate and optimal number of general practitioners (GPs) because of the high cost and long duration of medical education courses.
GPs are the primary providers of community health services and, as the referral agents, play an important role in satisfying the requirements of society (4). In the Iranian health system, due to the lack of a central registration system and the multiplicity of statistical centers, there is no precise data on the number of active GPs (5). It has resulted in concerns about the implementation of the Family Medicine Program (5, 6). Moreover, there is no mandatory referral system for optimal use of health services. Therefore, patients are more inclined to use the medical services of the specialists and sub-specialists for simple illnesses, especially given the low difference between patient pay for the GPs and the specialist services, which has enhanced their motivations. On the other hand, since the tariffs for specialist physicians are also fixed, they also prefer to visit patients with simple illnesses, as it takes less of their time (7). Therefore, when the GPs and specialists are compared, the demand for GPs has been less than the demand for specialist physicians.
In addition to the competition among physicians, the supply and demand of the GPs are influenced by some other various factors, which can lead to a shortage or surplus of GPs. The physician demand is generally influenced by national income per capita, out-of-pocket (OOP), age structure, and others (1, 3, 8, 9). Besides, the physician supply is also affected by national laws for increasing or maintaining the capacity of medical education, migration and Health brain drain, and the death and retirement of the physicians (3, 10, 11). Therefore, due to the need for careful planning in human health capital, especially the GPs as the first providers of the health services and referring factors having a direct impact on the health status of individuals, the impact of the imbalance in the GPs supply on the economy and health status of the community, and the need to adopt appropriate policies regarding the shortage and surplus of GPs have been investigated.  This study aims to estimate and predict the shortage and surplus of GPs in Iran and answer the following questions: What factors can influence the demand of general practitioners in Iran? Furthermore, are we facing a shortage or surplus of GPs in the future?
Numerous studies have been conducted on estimating the health workforce using different methods of need (12-16) and demand (8-10, 17, 18), some of which were cross-sectional, and some others were in time series. In recent studies, Scheffler and Arnold (3), Liu et al. (1), Dall et al. (19), and Chojnicki and Moullan (10) started to estimate and predict the physician while developing the conceptualized model for the physicians' demand by national statistics and time series. In line with these studies, the present study estimates the supply and demand functions of the GPs in Iran with time-series data at the national level and predicts the surplus or shortage of the GPs. Using Autoregressive Integrated Moving Average (20) and Vector Error Correction (VEC) econometric method to estimate and predict GPs supply and demand, perceptively, is the main difference between the present study and other studies.
The supply of health workforce includes all people who have the necessary expertise and skills to work in the health labor market (1). The supply of the health workforce is always affected by various factors, which can be classified into two categories of entry into the labor market (immigration, graduates, wages, and others) and exit from the labor market (retirement, migration, death, and others) (21). Salary is another important factor for entering the healthcare job market. The higher wages are, the more willing the specialist to be employed will be in the labor market (1, 3). High school students are also more encouraged to study in the health-related fields. As a result, the number of entries to medical universities and the supply of health-care workforce are increased (1). Also, many students are studying experimental science in Iran to enter the medical universities, compared with the mathematics and technical majors. Although a flexible wage is needed to balance the supply and demand of the health workforce (22), wages in the health sector are often not flexible because legislators initially set health workers' wages in most countries. Secondly, the health professional workforce has the great bargaining power to raise wages, and third, it is not easy to employ health workers in specific places for low-wage (3). Immigration and emigration of the health workforce have a huge impact on the supply of health workforce, especially health professionals, such as physicians (3, 10). African countries have the highest rates of physician migration due to the better benefits physicians receive when they migrate (23, 24).
The demand for the health workforce indicates a willingness to pay for health care by buyers (government and the private sector) (1). In other words, demand for the health workforce is a derived demand from health care and services. As the demand for health care and services grows, so does the demand for the health workforce. Therefore, the demand for health workforce is influenced by factors such as family income, financial ability of the government, conditions of demographic and epidemiologic and health coverage levels. (1). Increasing family income, government expenditures in the health sector, and aging are expected to increase demand for the health workforce (1, 8-10).
Also, increasing the level of health coverage will lead to an increase in the demand for health labor because people are more able to buy health services (3, 25). There are several ways to predict the physician demand, including the need, demand, population ratio, and service target-based approaches (12, 14-16, 26, 27). Each of these approaches has its advantages and disadvantages (16, 26), of which usually the combination method is the most suitable one (10, 16). The developed models by Liu et al. (1), Scheffler and Arnold (3), Chojnicki and Moullan (10) were used to predict GPs demand in this study, even though the quantity and quality of GPs' services, epidemiological factors, innovation, and new technologies have not been taken into account in those models (28). However, they are the perfect models for functionally and theoretically expansion.
Methods
The present study is an analytical and applied study conducted nationally using ARIMA and VEC models. The statistical population was Iran, and sampling was not done. The data were annual time series, which were extracted from various databases. GDP per capita, Out of Pocket Payments Ratio (OOP), Ratio of people over 65 to 14-65 years (AGING) were extracted from World Bank databases (24). Moreover, general practitioners per 1000 population (GPs) and the number of hospital beds per 1000 population (BED) were extracted from the yearbooks of the Statistical Center of Iran for the years 1991-2017. Due to the nature of the time series variables, stationary of the variables was necessary to avoid the spurious regression. Augmented Dickey-Fuller (ADF) unit root test for stationary of the variables was used. In this test, the null hypothesis indicates the unit root or non-stationary of the variables. Therefore, if the calculated value of the test was greater than its critical value, the null hypothesis based on the non-stationary variable was rejected. Next, the models were used to estimate the supply and demand of GPs, and how to estimate its coefficients were described in separate sections. The models and required tests were estimated in Eviews 10 software at a 0.05 significance level. Ethical approval was not required for this study because it is based on secondary data that are fully anonymous and publicly available.   
The ARIMA model for Supply of GPs
The past trend is generally used to predict the supply of the physician (1, 3, 9, 10). It is usually appropriate for these predictions to use univariate prediction models, including ARIMA, which are frequently used (29). The ARIMA method is one of the popular time series methods and a combination of three methods, including Auto Regression (AR), Moving Average (MA), and the integration, which are generally referred to as ARIMA (p, d, q) (26). P represents the order of the AR method and stands for the number of autoregressive terms, which associates the given variable with its previous values. Q is the order of the MA method and states that the desired variable is related to lagged error terms q times, and d denotes the integration order of the variable. The general form of the ARIMA model (p, d, q) is as follows:
               
 Error terms. If the order of integration y equals one, the y variable enters the equation as the first difference. The ARIMA method was used to predict GPs supply according to Box-Jenkins methodology in this study. Box-Jenkins proposed four steps for the ARIMA method (29). In the first step (identification), the ADF unit root test was used to determine the integrated degree (d) of the GPs. To determine p and q values, correlogram and partial correlogram were used. Therefore, at this stage, the ARIMA (5,1,1) was chosen. In the second step, the previous step model was estimated by the maximum likelihood method. In the third step, all the relevant models (25 models) were estimated to determine the best model, and ARIMA (5,1,1) model was selected because of the lowest degree of Akaike's Information Criterion (AIC). Finally, this model was used for forecasting within the sample period and outside the sample period in the fourth step.
The VEC model for the demand of GPs
In this study, the developed models by Liu et al. (1), Scheffler and Arnold (3), Chojnicki and Moullan (10) were used to predict GPs demand. At the macroeconomic level, most studies used Gross Domestic Product (GDP) and Gross National Income (GNI) to predict physician demand (8-10) because physician training courses are costly, and only countries with proper economic growth can invest in physician training. Demographic and health coverage measures were also entered into the model as factors influencing the physician demand (1, 3, 8). For the demographic measure, the ratio of the number of people over 65 to the number of people aged 14-65 years was used, and it is expected that by increasing this proportion, the demand for the health workforce increases. Many studies have also shown the direct influence of age structure and aging on physician demand (9, 10, 12, 14, 26). For the health coverage measure, the ratio of household out-of-pocket (OOP) payments to health expenditures was used. OOP payments for health care expenditures indicate the community's social support for the household health expenditures. Since the health insurance coverage, government subsidies, and other financial supports affect the OOP payments (1), the lower the level of health insurance coverage is, the greater the OOP payments will be. It is expected that using health services and consequently the demand of the GPs will decrease. Other factors that may affect GPs demand are the health resources used in the present study by the number of hospital beds per thousand populations, and it is expected that with the increase of this proportion, the GPs demand will increase as well. Additional structural factors, including training capacity, labor regulations, migration, and others, may also affect the labor market (1). However, they were not included in the model due to insufficient data. Therefore, the conceptual model of the present study was as follows:
GPs: General practitioners per 1000 population
LnGDP: Natural logarithm of GDP per capita based on purchasing power parity
OOP: Out of Pocket Payments Ratio for medical expenditure to total health expenditures
AGING: Ratio of people over 65 to 14-65 years
BED: The number of hospital beds per 1000 population
There are various methods for estimating and predicting the above model; however, since we expect GPs demand to be influenced by the current quantities of the independent model variables, they are also affected by their lagged values. The current GPs demand may also be affected by the demand for GPs in past years. Therefore, using a Vector Autoregressive (VAR) approach allowing for the indigenousness of all the variables in the model and the effect of different interruptions of the variables on one another seems good. On the other hand, since all model variables are time series and the possibility of a unit root and non-stationary that can lead to spurious regression and misleading results, it is necessary to use co-integration methods Vector Error Correction (VEC) models. Due to the co-integration of the model variables, the VEC model was used, the general structure of which for the model variables is as follows.

ECM: An Error Correction Component that represents a deviation from a long-run stochastic trend
D: first-order difference of the variables
P: Optimal lag
The Johansen method (30) was used to estimate the above equations and the long-run relationship between them. In this method, first, the optimal lags of the model should be selected based on the AIC, Schwartz Criterion (SC), Hanan-Quin (HQ), the final prediction error (FPE), and the LR test (LR), the optimal lags must have the least of these criteria. Secondly, the number of co-integration vectors is determined based on the maximum eigenvalue and Trace tests. In these two tests, the null hypothesis is the number of co-integration vectors, and if the value of these two tests is greater than the critical value, the null hypothesis will be rejected. After performing the above steps, the study models were estimated using VEC models (32, 33) in Eviews 10 software, and based on the results, the GPs demand was predicted.
Results
From 1990 to 2017, the GPs variable per 1000 people had a mean of 0.3022 ± 0.030 and moderate sinus growth; on average, its growth rate was 2.45%. Also, according to the Jarque-Bera statistic had a normal distribution. As a percentage of total health expenditures, the OOP payment variable had a slight downward, from 59.60% in 2000 to 38.79% in 2017, averaging 52.24 ± 7.07%. The number of hospital beds also had a slight upward trend per 1000 persons, which averaged 1.64 ± 0.13 over the period 1990–2017. This variable has multiplied since 2013. According to Table 1, ADF results showed that all the variables were non-stationary at the level, which was stationary at the first difference.

Projecting supply for GPs
In the first step, the results of the ADF test showed that the non-stationary GPs variable became stationary by a one-time difference. Thus, it is integrated of order 1, I (1), (Table 1). Since the variable of GPs supply was I (1), the first-order difference of GPs variable was used to investigate the correlogram. The correlogram and partial correlogram results showed that AR (5) and MA (1) might likely exist. Therefore, at this stage, the ARIMA (5,1,1) was chosen. In the second step, the previous step model was estimated by the maximum likelihood method. In the third step, all the relevant models were estimated to determine the best model (Table 2), and the ARIMA model (5,1,1) had the lowest AIC (-6.06), so it was selected. The Jarque-Bera statistic was 0.87, with a probability of 0.61, indicating that the error terms were normal. Also, Kurtosis and skewness of the error terms were normal. Unit root test for error terms also showed that the value of ADF statistic was -4.6021, which is (1% significance) higher than the critical statistic value (-3.72). Accordingly, the null hypothesis or non-stationary of error terms is rejected, and error terms were stationary; then, this model was selected as the best model used in the forecasting stage. The ARIMA model was first used for good prediction within the sample period (1991-2017), and then, the prediction results were matched with actual GPs supply results (Figure 1). It was observed that the prediction values are very close to the actual values, so this model was used to predict outside the sample period (2018-2030). The results of the supply projection from 2018 to 2030 were presented in Table 3. Also, the results of supply and demand projection for GPs were presented in Figure 2. The average GPs supply was estimated to be 0.2022 during the period 2018-2030.

Estimating and projecting demand for GPs
Models (3-6) were used to estimate and predict GPs demand. To avoid spurious regression results, the stationary of the variables was estimated. According to the results of Table 1, GPs, lnGDP, OOP, AGING, and BED are non-stationary at the level and stationary at the first-order difference, i.e., they are integrated with I (1), so there is a possibility of co-integrating vectors between them. Thus, Johansen co-integrating method (30) was used (32, 33). In this method, the optimal lags and the number of co-integrating vectors should be determined, and then the desired model is estimated using this information. According to AIC, SC, LR, FPE criteria, the optimal lags must have the least of these criteria.
Consequently, the first lag was determined as the optimal model lag (Table 4). Trace and maximum eigenvalues tests were performed to determine the number of co-integrating vectors. The results were presented in Table 5 and showed that the hypothesis of one co-integrating vector is not rejected at the 0.05 level. Therefore, there is at most one co-integration vector between the model variables. Due to a co-integration vector, the VEC models were used to estimate the models and investigate the relationship between the variables.

The long-term and short-term results of GPs demand were presented in Tables 6 and 7, respectively. For stability of the model coefficients, the inverse roots of the model were estimated, all of which were less than one and within a single circle. The error terms of the above model were tested for stationary by an ADF test. The statistic value was -4.36, which was greater than the critical value of -3.95 at the 5% level, so the error terms were stationary. The normality test of error terms showed that the chi-square values for Kurtosis and skewness were 4.05 and 2.36, respectively. Besides, the Jarque-Bera statistic value was 6.42, all indicating that the error terms were normal. The LM test of the residuals' Autocorrelation also showed no serial correlation at the model. Also, VEC residual Heteroskedasticity tests showed that the chi-square statistic value is 125.28, which is accepted the variance homogeneity with a probability of 0.36. Thus, the results of all diagnostic tests after estimating the model indicate the stability and reliability of the model results. Hence, this model was used to predict the number of GPs. Firstly, within the sample, forecasting was performed. Regarding the high consistency of the GPs prediction values and their actual value, this

model was used for outside sample prediction from 2018 to 2030, and the results were presented in Table 3.