A generative model for age and income distribution

Each individual in society experiences an evolution of their income during their lifetime. Macroscopically, this dynamic creates a statistical relationship between age and income for each society. In this study, we investigate income distribution and its relationship with age and identify a stable joint distribution function for age and income within the United Kingdom and the United States. We demonstrate a flexible calibration methodology using panel and population surveys and capture the characteristic differences between the UK and the US populations. The model here presented can be utilised for forecasting income and planning pensions.


Introduction
A universal element of societies is the emergence of hierarchical organisation structures within professions. People develop work experience through time and manage to obtain jobs of increasing responsibility and increasing level of income with time. Hence, it is a natural property of income distribution to be correlated with work experience and age; nevertheless, most income models do not study the relationship between income and age, and consequently between income distribution and demographic changes. This paper introduces a model of income, dependent on age-specific model parameters and random shocks. The model contributes to the understanding of the relationship between age and income and its dynamics.
Our aim is to compare the estimated parameters in the UK and the US age and income distribution to find out similar characteristics of age and income across states, as well as the contrasting differences. A simple age and income model is fundamental for the development of a sustainable pension system. The model focuses on the age and income relationship and further factors, such as occupational levels, are not considered. The model is estimated via panel survey data from the UK and population survey data from the USA. The data from panel surveys track the same individuals for the duration of the survey, and the population survey is repeated with different people each wave. The results reflect a clear income-age relationship in the UK and US, a clear structure of the joint distribution characterised by rapidly increasing income at younger ages, followed by income levels stabilising near mean income but spreading till retirement. At this point, the income decreases and concentrates around mean retirement income. The paper demonstrates a flexible methodology to estimate parameters from population surveys, as well as panel surveys. The paper provides a simple generative model to evolve age-income population for simulation and forecasting purposes, which can constitute the foundation for future studies of financially sustainable pension systems by providing a benchmark for capturing age and income relationship. The purpose is to have a baseline model simple enough for isolating age and income relationship of income dynamics. Such a model will serve for investigating the properties of a sustainable and balanced pension system. The mean and standard deviation statistics from the panel and population surveys on Fig. 6 and Fig. 1 from observed panel data and simulation results reflect a clear relationship between age and income. More complex models, which investigate additional factors, and profile heterogeneity of income dynamics are out of the scope of our work.
Previous research on income have been conducted, and the research focused on investigating and explaining wage dynamics. Champernowne explicitly introduces a first-order Markov process to model the time-evolution of wages [1]. Following Markov process path, the validity of the first-order Markov assumption is tested by Shorrocks [2]. Following research introduces a second-order Markov process, yet neither of these works links individual wage dynamics to time-evolution of the distribution of wages [3]. A different approach focused on poverty, which deals with modelling individual data using linear regression and transitioning to poverty (probit model) [4]. A more comprehensive model incorporating various factors is developed to estimate transition probability in wage quintiles conditioned on various regressors, including education, experience and age [5], furthermore study both intra-and inter-group inequality. The persistence of the low pay state and factors affecting the low pay probability are expressed with a generalized regression model. For modelling low income transitions the previous research use British Panel Data for the 90s, focus on the transition probability and state dependence for the poverty status [6] and define poverty transition equation as coarse-grained dynamics. Inequality and upward mobility between quintiles considering gender effects are investigated [7]. The previous models in literature either incorporate numerous external variables, distribution characteristics and functions, such as innovation constants or they are limiting their scope to the investigation of dependence on a single variable [8,9]. A more recent article by Guvenen investigated a model for which focal variables are the human capital consisting of education, work experience, and idiosyncratic shocks [10], following research modelled male income for studying the impact of labour income taxation policy on inequality [11] The referred life-cycle model's distribution characteristics of the pre-tax income arise from the differences in the individual's ability to learn new things and idiosyncratic shocks. Previous research tried to capture the income dynamics with Markov Models, linear autoregressive models, or by relying on econometric toolset such as covariance matrices. We investigate a generative model with an empirical distribution for sustainable age and income relationship in a population; we achieve this via an income evolution model with an age-dependent parameter, estimated from previous population and panel surveys.
The previous research [12] presents a two-class distribution, majority of population is described by the exponential function and small fraction of population of higher income individuals are described by power law [13]. The BHPS and IPUMS CPS data are topcoded and not suitable for studying the power law at the top, but the majority of the pop-ulation as reflected on Fig. 3 is consistent with empirically well-established exponential distribution of income [14].
Although there are models that incorporate indirectly the age as years of experience in job for studying income dynamics. There are no studies, to the best of our knowledge, studying the joint distribution of age and income in the scope of income evolution.
In contrast to previous research, our study introduces a dynamic model that describes the income-dependent only on age and previous income. This paper investigates the stationary property of the income distribution dependent on age. We provide a model in which the mean and variance of income given age are preserved at any time point.

Methods
We introduce a simple model which focuses on age and income relationship and differs from recent literature by not incorporating other variables such as occupational level, level of education and skill coefficients. The model is stationary, i.e. the mean and variance of income given age are preserved in time. The model is utilised to represent observed panel data for gaining empirical insights regarding age dependant, income dynamics and mobility. The calibrated model can be utilised as a simple generative model to evolve an age-income population for simulation purposes and it provides a theoretical background for studies focusing on ageing and pension income of the population. We initially assume the following model, by which μ(·) and σ (·) represent a function of age, income, and individual-specific additional parameters θ i or λ i , for the sake of generality. μ(·) is a function capturing mean income characteristics, and σ (·) captures the variational characteristics of the income. We consider the following individual income stochastic process for an economic agent i characterised at each time step t by a given age a i and income y i : The characterising insights on Fig. 1 from the panel data lead to the assumption that the probabilistic step at time t depends only on the age and income of the preceding step.

Defining income and age dynamics
Earnings of individual i at the time step t is denoted as Y i,t and its logarithm is y i,t . The parameters that describe the income process are: age-dependent persistence parameter q a , age-dependent mean μ a and age-dependent standard deviation σ a . The income shock process consists of independent random shock η i t which is normally distributed with mean zero and variance 1, and it is applied to σ a , the model can be defined as follows: Averaging income y, for individuals who are a years old, givesȳ a ,which denotes the average income for age group a across all individuals i and periods t. Assuming that the age-dependent income profiles are stationary, we can average incomes y i a,t across individuals and time to get the following equation: whereȳ a denotes the average income for age group a, taken across all individuals i and periods t. The following equation can find the estimator for μ a : The income data from different waves are inflation adjusted to isolate effects of economic growth.

Data
The British Household Panel Survey (BHPS) [15] from the UK and The Current Population Survey (IPUMS CPS) [16] from the USA are used for estimating the parameters of the model in Eq. (2) and comparing the results of simulated data and surveys. The BHPS is a Panel Survey conducted between 1991-2008. For our model we focus on labour income data, which captures wage, salary or self-employment income. To investigate population characteristics, we also incorporate other income sources and call it "Total Income", which additionally captures the transfers, pensions, grants, aids, state-benefits, dividends, capital income and rents. BHPS provides individuals specific longitudinal weights for ensuring the representativeness of the population. Two types of weights are provided with BHPS. The first wave is weighted for adjusting population marginals at the households and poststratified to the population age by sex marginals. Consecutive waves are re-weighted to take into account sample attrition, variables such as address change, household region, age, sex, race, employment status, income total and composition, educational qualifications [17]. Panel Survey is conducted via questionnaires with tracked individuals of the initial sample. Detailed explanation of the relevant variables from BHPS dataset can be foınd on the Appendix. There is an extension to the sample population in 1999. For the USA, IPUMS CPS is used, which is annually conducted with different samples each year. In contrast with BHPS, the Labour Income from IPUMS CPS does not include self-employed income, and the weights are cross-sectional. Income distribution, age distribution and income-dependent age distribution from the surveys are utilised for parameter estimation and further analysis. q a , μ a , σ a are the key parameters estimated according to the proposed model. Following investigation and interpretation of the estimated parameters, these parameters are used to simulate the population's income transitions. The simulation is initialised using the panel data from the Wave 1, and the income evolution function on Eq. (2) is applied transitively in an iterative approach to the data for simulating successive waves. The simulated data is plotted, interpreted and finally compared with the observed data. Figure 2 reflects the Population Pyramid in the UK, and how the shape evolved over the 18 years considered. The population pyramid of USA can be found in the Appendix. The UK population sample from BHPS has a relatively balanced population with a slight weight towards younger cohorts initially in 1991, which denotes Wave 1. The UK population gradually got older, and the population pyramid reflects mass's shift towards older generations, this shift happened gradually over the years. The US population from CPS reflects a young population in Wave 1 with a notable skew towards younger cohorts, after 17 years the US population loses this property towards younger cohorts and gets significantly older. Both the UK and US population get older and reflect a trend towards an ageing population, which will significantly impact the pension system.
The shape of the population pyramid and its evolution with time from the panel survey reflect an ageing community [18]. JDFs of Total and Labour Income in the UK and USA reflect that, there is a sharp increase between the ages 15-20, which can be interpreted as the beginning of the work-life, transitioning from part-time work to full-time work, and graduation from higher or vocational education. The most significant difference of the UK Total JDF in contrast to the UK Labour JDF is the tail section corresponding to the retired population, which denotes the significant percentage of individuals older than 55. The tail section is relatively concentrated, which can be explained by the state pension benefit levels and mandatory social security system. The US population reflects a surprisingly sparse older cohort for the Total Income data, and the most significant difference to the UK is the relatively lower income levels compared to the wage income. In the US population higher variance spread to a wider band, which might be caused by a non-standard retirement system not supported by strong state pension benefit and mandatory pension schemes during employment.
The comparison for the model simulation and observed data shows common characteristics, as the joint distribution of age and income in logarithmic scale is presented in Fig. 1: an initial sharp rise between the entrance to the graph on 16 years old, the amount of 16 years old includes pocket money, allowance and part-time or internship jobs. There is a steep increase in mean and variance between the initial income and income at the age of 20. The increase is sharper for the mean in comparison to the variance. The population's mass has similar characteristics with near 23k GBP annual income, and for ages between In compliance with the literature [12], the population is divided as two-classes, and the majority of the population covering low and middle income follows Boltzmann-Gibbs distribution. The observed and simulated populations from the UK and USA are reflecting exponential characteristics for low and mid-income individuals, log-linear PDF plots reflect similar PDF characteristics on Fig. 3.
In the following sections we will focus on labour income and employed labour population. Total Income covers all of the income streams including transfer income such as pocket-money, labour income, capital income, and pension income; these different streams might be governed by varying dynamics non-uniform across the type of income; so we decided to focus on labour income, which involves the broadest section of the population; with most significant impact. The only other primary source of income in terms of gross value is the capital income, which might be significantly affected by other factors such as inter-generational shifts, market conditions and global financial state. In order to focus on labour income dynamics, the other income sources are left out of our modeling.

Data processing and calibration
BHPS provides a vast amount of socio-economical data for each individual and household participating in the study. The columns of income, age data, the individual's statistical weight representative of the British population and overall survey-with the individual's intra-wave unique identifier mostly suffice for this paper's scope. PID, Wxage12, wFIYR, Wxrwght fields of BHPS are used for each wave.
The Income variable xfiyr is each individual's annual income, including labour income, benefits, pensions, transfer income, and investment income. Participants were asked according to annual income in the reference year from September in the year prior to the interview until September in which interviewing begins [17]. The income figures are adjusted for inflation, as part of pre-processing. During the dataset preparation, a floor wage is determined to exclude in labour income, which denotes to excluding part-time and short-employment income. The income data is inflation-adjusted and transformed into log-domain.
IPUMS harmonizes the CPS and provides IPSUM CPS micro-data. The IPSUM CPS includes a large spectrum of topics such as demographics and employment, as well as supplemental studies such as the Annual Social and Economic Supplement (ASEC). Each individual can be identified by "CPSIDP", "INCTOT" and "INCWAGE" correspond to the total income and wage income, and "ASECWT" denotes the weights derived from ASEC Supplement. The data set is topcoded, and specific codes are used for labelling missing and incorrect data. The ages over 80, 90 and 99 are top-coded till 2004, and after 2004, the top-coding bins are determined as ages 80, 85, 90 and 99 by the panel data collectors [16]. Although this dataset contains high-income individuals, there is top-coding applied, so individuals with very high income are not included.

Fitting distributions
Estimating the income evolution function parameters is the most critical part of the research, and the decision depends on various factors such as the type of data, bias, and assumptions. Various techniques are investigated, leading to different results, with each having unique strengths and weaknesses.
The first method investigated is Generalized Method of Moments (GMM), which presumes that the first three moments of the income evolution functions provide the necessary information for approximating the underlying generative process. The equations of the first three moments of the income evolution equation can be solved for the parameters Figure 4 q a , σ a and μ a Plots for UK Labour Income, the parameters are estimated with two different methodologies of Generalized Method of Moments and Least Squares Method. The LSM utilises longitudinal data, which provides capability to estimate substantially high q values, which represent the persistence level of an individual's income, GMM capture much lower q values, due to its incapability of utilising longitudinal data q a , σ a , μ a . Both of the BHPS from the UK and the IPUMS CPS from the USA can be used for estimation with Generalized Method of Moments (GMM) with first three moments as reflected on Fig. 4 and on the Appendix.
The second method utilises the micro-data from the longitudinal surveys, which tracks the individual for consecutive years. The parameters are approximated to fit the income evolution function using least squares minimisation for the individuals participating in the studies for consecutive years. The BHPS from the UK is a suitable micro-data consisting of a panel survey, and the survey tracks income of the same individuals over the years.

Estimation for generalized method of moments (GMM)
The three moments of the age income evolution function are utilised to find a polynomial; afterwards, the equation is solved for q a , σ a and μ a ; at this point, a observed solution for parameters is found, but the relationship captures only the dynamics of the first three moments. Calculations can be found in the Appendix. The statistical variables such asȳ a are found for each wave and than averaged across waves for finding a one set of stationary variables, which can be used to estimate q a , σ a , μ a . The details and derivations for the GMM estimation technique can be found in the Appendix.

Estimation of least squares for micro data
The Least Square Method requires that an individual's income for two consecutive years be existent in the dataset, this restriction is fulfilled by the BHPS, a panel survey, but the CPS IPUMS population survey does not satisfy this condition. The income data from two consecutive years per agent is used to estimate age-specific parameters, which characterise the income evolution function at Eq. (2). LSM tries to estimate parameters by fitting the data to the income evolution function.

The generative model
The model can also be used for simulation and forecast, tracking income trajectories of the individuals, providing a bench table for observing the stylized facts and complex properties of the income dynamics. Following the estimation of model parameters, the model is bootstrapped with data from Wave 1 for initialising the simulation. Each individual from Wave 1 is initialised as an agent in our model. According to Age Income Evolution Dynamics Eq. (2), the income of an agent is transitively updated at each consecutive wave update. η it provides the random feed, which introduces variability for the income evolution of the agents. At each wave update, a new generation of agents consisting of 25 years old individuals from the initial wave are injected. Following each wave, distributions corresponding to the state of the simulated population are calculated. A full calibration of the model is shown in the Appendix.

Parameters
The optimized performance of these three methods are compared and discussed in the following sections. No boundaries are explicitly imposed by LSM estimation. The μ, σ and q variables are independent of each other, but the estimation process or data itself can introduce a slight dependence. The GMM estimation technique results in minimal q values near 0, so the estimated parameters approximately resemble an auto-regressive model. However, despite near 0 negligible q values, the q plot has a distinct shape with an increasing trend with a small decrease between 25 and 30, has very different characteristics depending on the estimation method. The GMM estimation method results in minimal q and the μ reflects the characteristics ofȳ, which is in compliance with this estimation method's nature. The μ value increases at first and then plateaus and slightly decreases near retirement. On the contrary LSM estimation mainly characterises the income with an increasing q parameter, so the μ parameter has limited effect and reflects a

GMM
GMM estimation technique approximates the μ a values to be consistently around 10 and the q a values are around 0 with an initial sine-like wave followed by a steady increase. The σ a values are around 0.8 and have a positive trend. q a values display a positive trend as well. Theȳ a and std(y) a plots of the simulation is similar to the observed data, but the standard deviation plot is particularly noisy. The JDF of the simulation on Fig. 5 is sparse, consistent; but not highly concentrated around mean. Both of these methods depends on assumptions about the dynamics of the income evolution function. The GMM method assumes that the first three moments of the equation are sufficient for estimating the parameters because they provide a solvable system. However individual characteristics in an age group such as different income levels and clusters within are lost during the moment estimation.

LSM by individual transitions
In order to use LSM for approximating the parameters, one needs the individual income transitions in consecutive years, thus identifying the same individual in consecutive co-horts is necessary and the panel studies such as BHPS satisfy this condition. The agedependent income evolution function is fitted with individual income transitions of consecutive years, which results in consistent parameter plots and theȳ a plot of the simulation reflect similar shape with the observed data on Fig. 6. The JDF of the simulation on Fig. 7 is able to reflect the dispersion among various clusters better because unlike the other methods heavily depending on the statistics such as mean and standard deviation of the entire age group, the LSM utilises individual-level microdata.
The 95% confidence interval with 2000 bootstrap samples for the estimated parameters from UK microdata by LSM can be found on Fig. 8. It is evident from the plots ofȳ a and std(y a ) for the observed and simulated data that the model can capture the characteristics of the income conditional on age distribution. A close investigation of simulations on models calibrated with UK Labour Income Data suggests that the GMM is most successful for reflecting the outcomes with similar mean and standard deviation characteristics of all waves after simulation with 18 waves that were simulated with the parameters q a , σ a , μ a estimated by the GMM. But LSM reflects the individual trajectories, and JDF more accurately. The results showing the performance of GMM method can be found on the Appendix.
A general analysis of the comparison of joint distribution of age and inflation-adjusted income results in the following plots for weighted observed data and simulated data in Fig. 5: JDF of the simulated UK Labour data is in parallel with the expectations for GMM Estimation method, consistent and stable, resembling a similar shape but not concentrated for the heat regions with intense concentrations on Fig. 6. The main difference between the observed and simulated JDFs is the concentration of the mass of the population between 23 and 50.

Wave-specific analysis
The population from wave-1 is used for bootstrapping the simulation and the weights of the individuals are not incorporated to the simulation, because the income evolution Eq. (2) is the focus of this paper, and the main purpose is not the perfect representativeness of the initial wave. The new agent injection on 1999 by panel survey is reason of difference in the UK simulation and observed JDF plots. Although the simulation's initial state is bootstrapped as the unweighted dataset, starting from the second wave, the JDF of the simulated population resembles the characteristics of the JDF from the panel survey with the weighted population, which reflects that the model is successfully capturing income evolution dynamics.

A simple pension system
A financially sustainable pension system can be characterised by the balance between inflow and outflow of funds. The specifics and stability of pension system is out of the scope of this paper, and needs case specific detailed modelling. For a general demonstration, we assume simple inflow and outflow dynamics(Eq. (5) and Eq. (6)), which are derived to represent statistical properties of the savings and consumption. Figure 9 reflects the imbalance between inflow and outflow, which results in a deficit.
Pension is assumed to be annually £16,368, in light of the median net income before housing costs for all pensioners from DWP Pensioners Income Series in 2008/2009 [19]. Constant alpha for pension saving rate is selected to be 0.0775 to 0.2, in light of OECD Pension Report statistics [20].
The amounts are adjusted for inflation and reflect the 2009 levels. The inflow and outflow plots from our simplified generalisation of the pension system reflect a deficit.

Discussion
The income evolution Eq. (2) of the proposed model consists of the parameters q a , μ a , σ a : the persistence coefficient for the respective age group q a , determines the rate of persistence at a given age.
Age-dependent mean income parameter μ a expresses the expected age-specific income evolution mean for the next income and behaves such that if the mean parameter is high the persistence q a is most likely to be lower. If the mean parameter μ a is lower, the persistence parameter is higher which signals a potential widening of the income gap for the population.
σ a captures the variability of the individuals according to conditional distribution and incorporates randomness of the shocks.
The social safety nets, basic pension incomes and the Defined Benefit Pension plans are financed via the working population; the ever-growing unbalance towards ageing cohorts needs careful forecasting and planning. The demographic shift will impact the economy's functioning in general, introducing a heavy burden to welfare states financing the health  and pension of the retired population, which will reflect society as taxes and benefit cuts. The best course of action is forecasting the changes and planning in advance for the future.

Interpretation
The q a persistence estimated by GMM reflects that UK population reflects an initially high q a value in youth, followed by relative decrease, and then a consistent increase. The q values estimated by GMM fluctuate around 0 and minimal. The income persistence variable of individuals is not captured by GMM, which does not utilise panel survey's tracked individual income micro-data each year.
LSM calibration methodology uses longitudinal information regarding evolution of an individual's income, this property of LSM makes it suitable to provide higher q values, which captures the persistence. The GMM calibration methodology does not utilise longitudinal data, which makes it applicable with survey data but results in lower q values. This tells us that we are still able to reproduce consistently the data year by year almost disregarding the past year data. Our analysis suggests that LSM is a more reliable calibration model as using longitudinal information appears to be crucial to capture the income evolution dynamics. Furthermore the model has power to capture some heterogeneity, if the parameters are fit with LSM; and the income evolution function can preserve the bootstrapped wave heterogeneity due to persistence parameter q and the randomness injected with σ can provide outlier behaviour after multiple waves of income evolution function updates.
The LSM results in a consistently increasing q a value by the UK model, with a significant jump between ages 25-30, which corresponds with a μ a plot consistently decreasing with a significantly sharper decrease between ages 25-30. q a and μ a corresponding each other in an inverse proportion, especially by significant changes, especially by LSM. There are various examples of parameter effects that can be observed from BHPS dataset.
One example of q a effect is the upward mobility of age-group between 25 and 30, which is reflected by the increasing q a values and sharp increase observed on the jointdistribution plot. This effect can be due to finishing higher education and internships, in addition few years of experience, which results in a widening of income scissors. This change in mobility is healthy for the economy and does not represent a negative effect. One assumption should be researched further; if either this initial difference in mobility might limit of people with lower income for upward mobility.
An example of the σ a mobility is the age group of 30-35, which is reflected by a locally sharp increase of σ a values. Such mobility reflects a bidirectional movement of income for individuals, and such a variation might arise from the short-time employment, interruption of employment for education, temporary jobs and most importantly this mobility might be caused by the initial differentiation according to the education of individuals such as higher education or vocational education. This window represents an increase in the variation of the income.
In general, the shape of the distribution can be explained in three periods; the first period is the introduction to employment and teenagers, which represents income from parttime and temporary jobs at the beginning and start of full-time employment it sharply increases on Fig. 1.
The age group of 25-55 denotes the main productive era of the economic life, and the income reflects a high dispersion. All of the factors and random shocks act together and result in dispersed but a consistent distribution. Mobility wise this era provides opportunities for upward mobility and possesses downward mobility risks. At the end of this period, income tends to decrease slightly, which reflects a decrease in productivity. Another limiting factor is the minimum wage and state benefits, which introduces a lower bound envelope for the mass. Income sources and affecting factors of individuals in this era vary greatly, which results in the widest dispersion in the entire life-span. Some of the factors are education, social strata, adaptability to innovation, total-work hours per week, experience and expertness, seniority of the jobs and ageism. The third and final era represents the exit from the workforce and retirement, and temporary or part-time jobs for low-income old individuals. The income decreases gradually as the number of individuals exiting workforce increases with time, the income stabilises, and variation decreases significantly. Income in this era is relatively low, and the source is usually pension benefits, state support or temporary jobs. This model's outcomes can be used for various purposes; the most apparent fields for drawing consequences and planning are the works on inequality and mobility depending on age. Characteristics of workforce entrance, work-efficiency of individuals per age, the structure of the society, pension system, income stability, and the taxation system are the most obvious fields.
In the paper, two main estimation techniques are investigated, and the corresponding results from the simulated waves are presented. The first estimation method investigated is GMM Estimation. The income regions appear smoothed and spread. The second estimation method investigated is LSM, it utilises the microdata and is suitable for capturing an agent's income evolution. The JDFs from the simulated waves have the most similar mean characteristics to the observed data.
The LSM evidently performs better by utilizing longitudinal microdata; the GMM estimation method can be applied to both population and panel surveys, provides feasible distributions but with unrealistic modeling of an agent's individual income trajectory.

Conclusions
We demonstrated (1) a clear income-age relationship, which is reflected by the data from BHPS and IPSUM CPS, as well as simulations. (2) a clear structure of the joint age-income distribution in both the UK and USA. (3) a flexible methodology to estimate parameters from population surveys, as well as panel surveys. (4) a simple generative model to evolve the age-income population with real constraints for evaluating general policy scenarios, that is agnostic about occupation levels.
The model can be interpreted as delivering a premise that the information of an individual's experience and education can be encapsulated by income. Although in early career, the income dynamics are governed by the initial difference at the level of education and profession; the main dynamics governing income transitioning can be reduced to the relationship between income and age, which collectively encapsulate education and experience. These premises can be leveraged for developing simplified models for evaluating mobility, inequality, welfare state, and pensions.
The proposed model focuses on the evolution of age and income population and the paper successfully demonstrates a simple model that can be calibrated for age and income that can be used as a backbone for forecasting income and planning pensions. Understanding the dynamics and having the ability to forecast the age and income population is the key to the design of financially sustainable pension systems.
There are different dimensions for the future work: one of the dimensions is injecting random shocks to the distributions itself, which can be in the form of new population injection or withdrawal, as well as tuning the η it with various means for simulating a global or regional shock, such as pandemics or mass migration. Stress-testing the age and income distribution for different labour market scenarios could lead to relevant policy implications. The second dimension for future work is modifying the simulation system to estimate parameters on the fly, and provide a more adaptive and granular version of the simulation system. The third dimension for future work is incorporating data encompassing more years and more countries and with a higher resolution in time to investigate the role of multiple economic factors for short, medium and long time horizons. A.1 Variables of BHPS dataset

A.2 Model calibration
We can define the mean and standard deviation of income at a given age a as following: y i a,t =ȳ a .
The standard deviation and mean has the following relation with the squared average of incomes: η it has characteristics of the standard normal distribution: Squaring both sides of income evolution equation (2) results in following distribution: Equation (9) can be formalized as: Placing Eq. (14) for a + 1 and Eq. (13) results in following equation: (ȳ a+1 ) 2 + std(y a+1 ) 2 = q a y i a,t + μ a + σ a η i t 2 .
Expanding the right side of the equation results in: (ȳ a+1 ) 2 + std(y a+1 ) 2 = q a y i a,t (ȳ a+1 ) 2 + std(y a+1 ) 2 = q a y i a,t 2 + μ a + σ a η i t 2 + 2 q a y i a,t μ a + σ a η i t (17) = q a y i a,t 2 + (μ a ) 2 + σ a η i t 2 + 2(μ a σ a η it ) + 2 q a y i a,t μ a + q a y i a,t σ a η i t .
This equation can be solved for q a corresponding each age group. Cardano solution for cubic equations guarantees single real root to exist, the other two complex roots that Cardano solution provides are not used. Both of the σ a = std(y a ) and GMM estimation techniques can use the following equations for determining the μ a and σ a : For (q a ) 1 and (q a ) 2 according to Eq. (4): The σ 2 a can also be expressed in terms of q a , using Eq. (4)): σ 2 a = ( a+1 ) 2q 2 a ( a ) 2 -ȳ a+1q a μ 2 2 -2q a (ȳ a+1q aȳa )ȳ a .
A.4 Supplementary plots of the USA Figure 10 Population Pyramid for the USA Income between Ages 15-100 reflecting changes of age distribution in 18 years, which reflects mostly an aging population