Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
82,221 characters · 14 sections · 29 citation commands
Flexible Bayesian Models for Time-Varying Income Distributions
\affil{School of Mathematics and Physics, University of Wollongong, Wollongong, New South Wales, Australia}
Keywords: Generalised Beta 2 Distribution; Random Walk Model; Markov Chain Monte Carlo; Horseshoe Shrinkage Priors; Posterior Probability of Stochastic Dominance; Aboriginal Population Subgroup.
The estimation of income distributions plays a central role in the measurement of inequality and poverty and, more broadly, in welfare comparisons across time and across populations. For useful overviews of the extensive literature on income distribution modelling, including alternative specifications, their properties, and estimation methods, see the monograph by KleiberKotz2003, the edited volume by Chotikapanich2008, and the articles by BandourianMcDonaldTurley2003 and McDonaldXu1995. The present paper focuses on inference for income distributions of population subgroups for which only a small number of observations are available in each year.
The availability of detailed survey data has transformed empirical research on inequality and poverty. These data provide a rich basis for tracking distributional change, evaluating policy reforms, and quantifying the persistence of economic disadvantage. Prominent examples include the Household Income and Labour Dynamics in Australia (HILDA) survey WatsonWooden2012HILDA and the Panel Study of Income Dynamics (PSID) in the United States McGonagle2012PSID. In Australia, HILDA has become a key resource for analysing how welfare, inequality and poverty evolve for different demographic and labour-market groups. As such datasets continue to mature and expand, they increasingly support inherently dynamic questions: whether inequality and poverty are rising or falling, and how future distributions might evolve under plausible trajectories.
Policy discussions typically focus on inequality measures, poverty indices, and welfare comparisons that depend on the full income distribution. These include widely used measures of inequality, such as the Gini coefficient and generalised entropy measures, including the Theil indices Cowell2011, as well as poverty measures such as the headcount ratio, poverty gap, and related indices FosterGreerThorbecke1984. In addition, distributional comparisons based on Lorenz and stochastic dominance provide partial orderings that are often more informative than any single summary measure barrett2003consistent,barrett2014consistent.
A common strategy in applied work is to estimate a parametric income distribution separately for each year, treating each year as an independent cross-section, for example Dagum Dagum1977, Singh-Maddala SinghMaddala1976, and Generalised Beta 2 (GB2) distributions McDonaldXu1995. Under this approach, the distributional parameters are estimated using only within-year observations. The posterior densities of inequality and poverty measures, together with the posterior probabilities of dominance, are then computed as functions of parameter draws obtained from Markov chain Monte Carlo (MCMC) algorithms gunawan2021posterior. While simple, this year-by-year estimation ignores an important empirical feature: income distributions typically evolve smoothly over time, interrupted by occasional shifts associated with macroeconomic shocks, labour-market changes, or policy reforms. When adjacent years are in fact similar, estimating each year independently can lead to noisy parameter paths and unstable year-to-year movements in derived welfare, inequality and poverty measures, particularly for small sample sizes. As a result, the independent cross-sectional approach may produce spurious volatility in estimates of inequality and poverty, and may also yield misleading posterior probabilities of dominance.
This paper makes several contributions. First, this paper develops a flexible dynamic statistical model for income distributions that explicitly links information across time. The key idea is to treat the distributional parameters as latent, time-varying states that evolve according to a stochastic process. Second, it shows that dynamic modelling improves inference for income distributions, especially for population subgroups with small annual sample sizes. By linking adjacent years, the proposed approach yields smoother and more stable estimates of income distribution parameters and substantially more precise inference for derived welfare summaries. Third, the paper shows that the benefits are not limited to narrower credible intervals. Since poverty, inequality, and dominance comparisons are sensitive to noise in the fitted distributions, independent year-by-year models can lead to substantively different and potentially misleading conclusions about welfare, poverty, and inequality changes over time. In contrast, the proposed dynamic models provide posterior inference that is more stable, and better aligned with the underlying temporal structure in the data.
We illustrate the approach using both a simulation study and an application to income data from HILDA, with particular emphasis on the Aboriginal population subgroup and the residents of the Australian Capital Territory (ACT) subgroup. The simulation results show that, relative to independent year-by-year models, the proposed dynamic specifications recover the underlying time-varying income distribution more accurately and provide tighter uncertainty quantification for poverty and inequality measures. In the empirical application, the dynamic models yield smoother parameter trajectories and more stable welfare summaries, and they can change posterior inferences about Lorenz and stochastic dominance over time.
The remainder of the paper is organised as follows. Section (\oldref{sec:model}) presents the proposed dynamic income models and discusses the prior distributions of the model parameters. Section (\oldref{sec:posteriorinference}) discusses Bayesian inference for the proposed models using Markov chain Monte Carlo. Section (\oldref{sec:inequalitymeasures}) discusses poverty and inequality measures. Section (\oldref{sec:stochasticdominance}) discusses Lorenz and stochastic dominance conditions. Section (\oldref{sec:simulationstudy}) discusses a simulation study. Section (\oldref{sec:empiricalapplications}) discusses the real data application using HILDA survey. Section (\oldref{sec:conclusions}) concludes. The paper also has an online supplement with additional technical details and examples.
In this section, we develop dynamic parametric models for income distributions observed over time. Rather than estimating each year's income distribution independently, which can lead to noisy parameter trajectories and unstable estimates of inequality and poverty measures, we treat each annual observed income as a realisation from an underlying distribution that evolves over time. The key idea is to allow the parameters of a flexible income distribution to vary by year, while borrowing strength across adjacent years to regularise estimation and improve robustness. This is particularly important for subgroup analysis, where some annual samples are small and year-by-year estimation is highly variable.
Let $\mathbf y_t=(y_{t,1},\ldots,y_{t,n_t})^\top$ denote the vector of observed incomes in year $t$, where $n_t$ is the number of individuals or households observed in that year. Collecting all years, $\mathbf y = \big(\mathbf y_1^\top,\ldots,\mathbf y_T^\top\big)^\top.$ We model the evolution of the income distribution over time through a set of time-varying latent state parameters.
This section discusses the proposed dynamic models for income distributions. Let $p_{Y_t}(\cdot \mid \boldsymbol\phi_t)$ denote the density of a parametric income distribution at time $t$ with the parameter vector $\boldsymbol\phi_t$. Examples include the Dagum, Singh--Maddala, and GB2 distributions discussed in Section (\oldref{sec:incomedistributions}) of the online supplement. Conditional on $\boldsymbol\phi_t$, the incomes observed in year $t$ are assumed independent and identically distributed: \[ y_{t,i}\mid \boldsymbol\phi_t \stackrel{\mathrm{iid}}{\sim} p_{Y_t}(\,\cdot \mid \boldsymbol\phi_t), \qquad i=1,\ldots,n_t,\quad t=1,\ldots,T. \] Hence the likelihood factorises as \[ p(\mathbf y\mid \boldsymbol\phi_{1:T}) = \prod_{t=1}^T \prod_{i=1}^{n_t} p_{Y_t}(y_{t,i}\mid \boldsymbol\phi_t), \] where $\boldsymbol\phi_{1:T}=(\boldsymbol\phi_1^{\top},\ldots,\boldsymbol\phi_T^{\top})^{\top}$. The income distribution parameters are often constrained; for example, scale and shape parameters are typically positive. To work on an unconstrained space, we introduce a latent state vector $\boldsymbol\theta_t = (\theta_{1,t},\ldots,\theta_{d,t})^\top \in \mathbb R^d$ and a smooth one-to-one transformation $g$, such that \[ \boldsymbol\phi_t = g(\boldsymbol\theta_t), \qquad t=1,\ldots,T. \] For instance, positivity constraints can be enforced using exponential transformations. Under this reparameterisation, the observation model becomes \[ y_{t,i}\mid \boldsymbol\theta_t \stackrel{\mathrm{iid}}{\sim} p_{Y_t}\!\left(\,\cdot \mid g(\boldsymbol\theta_t)\right), \qquad i=1,\ldots,n_t,\quad t=1,\ldots,T. \]
To model gradual changes in the income distribution over time, we place a dynamic prior on the latent states $\{\boldsymbol\theta_t\}_{t=1}^T$. A natural starting point is a random walk (RW) model. Assume that ${\phi_{1,1},\ldots,\phi_{d,1}>0}$ and define ${\boldsymbol{\theta}_1 = \bigl(\log(\phi_{1,1}),\ldots,\log(\phi_{d,1})\bigr)^\top,} $ with ${ \phi_{1,1},\ldots,\phi_{d,1} \stackrel{\mathrm{ind}}{\sim} \mathrm{Half\text{-}Cauchy}(0,1). }$ Under this log-transformation, the density of the initial state $\boldsymbol{\theta}_1$ is given by
where $\phi_{k,1}=\exp(\theta_{k,1})$, for $k=1,\ldots,d$. For $t=2,\ldots,T$, the state evolution is given by \[ \boldsymbol\theta_t \mid \boldsymbol\theta_{t-1},\boldsymbol\sigma \sim \mathcal N\!\left(\boldsymbol\theta_{t-1},\,\mathbf Q(\boldsymbol\sigma)\right), \] where $\mathbf Q(\boldsymbol\sigma)=\mathrm{diag}(\sigma_1^2,\ldots,\sigma_d^2)$ and $\boldsymbol\sigma=(\sigma_1,\ldots,\sigma_d)^\top$ controls the magnitude of year-to-year change in each latent state. Equivalently, componentwise, \[ \theta_{k,t} = \theta_{k,t-1} + \sigma_k \varepsilon_{k,t}, \qquad \varepsilon_{k,t}\stackrel{\mathrm{iid}}{\sim}\mathcal N(0,1), \qquad k=1,\ldots,d,\quad t=2,\ldots,T. \] This formulation allows different parameters of the income distribution to evolve at different rates. The innovation scales determine how strongly the model borrows information across adjacent years. Small values imply smoother trajectories, whereas larger values allow more pronounced temporal variation. We use a weakly informative half-Cauchy prior: ${ \sigma_k \sim \mathrm{Half\mbox{-}Cauchy}(0,1), }$ for $k=1,...,d$. Because the half-Cauchy prior is concentrated near zero but has heavy tails, it encourages small value of $\sigma_k$, while still allowing for large value of $\sigma_k$ when warranted by the data. Combining the observation and state equations yields the joint posterior density \[ p(\boldsymbol\theta_{1:T},\boldsymbol\sigma\mid \mathbf y) \propto \left[ \prod_{t=1}^T \prod_{i=1}^{n_t} p_{Y_t}\!\left(y_{t,i}\mid g(\boldsymbol\theta_t)\right) \right] p(\boldsymbol\theta_1) \left[ \prod_{t=2}^T p(\boldsymbol\theta_t\mid \boldsymbol\theta_{t-1},\boldsymbol\sigma) \right] p(\boldsymbol\sigma), \] where $\boldsymbol\theta_{1:T}=(\boldsymbol\theta_1^\top,\ldots,\boldsymbol\theta_T^\top)^\top$. This makes clear how the dynamic model differs from year-by-year estimation: inference for $\boldsymbol\theta_t$ is informed not only by the data observed in year $t$, but also indirectly by neighbouring years through the latent trajectory. As a result, the model yields smoother and typically more stable estimates of time-varying income distributions and of the inequality and poverty measures derived from them.
The random walk model can be extended by considering a horseshoe shrinkage prior carvalho2010horseshoe on the innovation scales. This yields a {random walk income model with horseshoe shrinkage priors} (RW-HS). The main idea is to replace the innovation scales by a global-local shrinkage structure: a global parameter controls the overall amount of temporal smoothing, while local parameters allow individual components of the latent state to deviate from this common level of smoothness when supported by the data.
Under the horseshoe specification, the latent states evolve according to (ref) for $t=1$ and, for $t=2,\ldots,T$, \[ \boldsymbol\theta_t \mid \boldsymbol\theta_{t-1},\tau^2,\boldsymbol\Lambda \sim \mathcal N\!\left(\boldsymbol\theta_{t-1},\,\tau^2\mathbf\Lambda\right), \] where $\mathbf\Lambda=\mathrm{diag}(\lambda_1^2,\ldots,\lambda_d^2). $ Equivalently, componentwise, \[ \theta_{k,t} = \theta_{k,t-1} + \tau\lambda_k \varepsilon_{k,t}, \qquad \varepsilon_{k,t}\stackrel{\mathrm{iid}}{\sim}\mathcal N(0,1), \qquad k=1,\ldots,d,\quad t=2,\ldots,T. \] Here, $\tau>0$ is a global shrinkage parameter and $\lambda_k>0$ is a local shrinkage parameter for the $k$th latent state. When $\tau$ is small, the latent trajectories are strongly smoothed across time. When a particular component requires more flexibility, a large value of $\lambda_k$ allows that component to have a larger jump. We assign half-Cauchy priors to both the global and local scales: ${\lambda_k \sim \mathrm{Half\mbox{-}Cauchy}(0,1), k=1,\ldots,d,} $ and $ \tau \sim \mathrm{Half\mbox{-}Cauchy}(0,1). $
This section describes the Metropolis-within-Gibbs (MwG) samplers RobertCasella2004 for the proposed dynamic models for income distributions.
We first consider posterior inference for the random walk model for income distributions. Algorithm (\oldref{alg:mwg_rw_income_model}) gives the full MwG sampler. Each MwG iteration alternates between two steps: (i) updating the latent states $\boldsymbol\theta_{1:T}$ one time point at a time using random-walk Metropolis-within-Gibbs updates, and (ii) updating the innovation variances $\boldsymbol\sigma^2$ using Metropolis-within-Gibbs steps with inverse-gamma proposals.
We define the time-$t$ log-likelihood contribution as $ \ell_t(\boldsymbol\theta_t) = \sum_{i=1}^{n_t} \log p_{Y_t}\!\left(y_{t,i}\mid g(\boldsymbol\theta_t)\right). $ At each iteration, we update $\boldsymbol\theta_t$ conditionally on its neighbours using the Gaussian proposal $ \boldsymbol\theta_t^\star \sim \mathcal N\!\left(\boldsymbol\theta_t,\kappa_t\boldsymbol\Sigma_t\right), $ where $\boldsymbol\Sigma_t$ is a positive definite proposal covariance matrix and $\kappa_t>0$ is a tuning constant. In practice, $\boldsymbol\Sigma_t$ is obtained from the empirical covariance matrix of previous MCMC draws, and $\kappa_t$ is tuned to target an acceptance probability around $0.2$ using the algorithm proposed by garthwaite2016adaptive.
Because the latent process is first-order Markov, $\boldsymbol\theta_t$ interacts with the rest of the trajectory only through its immediate neighbours $\boldsymbol\theta_{t-1}$ and $\boldsymbol\theta_{t+1}$. Hence, when proposing $\boldsymbol\theta_t^\star$, the acceptance probability depends only on $\ell_t(\boldsymbol\theta_t)$ and the transition densities linking $\boldsymbol\theta_t$ to adjacent states. Define \[ \log \pi_t^{\mathrm{RW}}(\boldsymbol\theta_t)=
\] Since the proposal is symmetric, the acceptance probability is ${ \alpha_t = \min\Bigl\{1, \exp\bigl( \log \pi_t^{\mathrm{RW}}(\boldsymbol\theta_t^\star) - \log \pi_t^{\mathrm{RW}}(\boldsymbol\theta_t) \bigr) \Bigr\}. }$ We then set $\boldsymbol\theta_t= \boldsymbol\theta_t^\star$ with probability $\alpha_t$. Detailed updates for the innovative variances are given in Section (\oldref{sec:additionaldetailrandomwalkincome}) of the online supplement.
We now extend the MwG sampler to the random walk income model with horseshoe shrinkage priors carvalho2010horseshoe. Under this specification, the state innovations satisfy \[ \Delta\theta_{k,t} = \theta_{k,t}-\theta_{k,t-1}, \qquad t=2,\ldots,T,\quad k=1,\ldots,d, \] with \[ \Delta\theta_{k,t}\mid \tau^2,\lambda_k^2 \stackrel{\mathrm{ind}}{\sim} \mathcal N(0,\tau^2\lambda_k^2). \]
A convenient augmented representation of horseshoe priors is obtained through inverse-gamma mixtures makalic2015simple. Specifically, \[ \lambda_k^2 \mid \nu_k \sim \mathrm{IG}\!\left(\frac{1}{2},\frac{1}{\nu_k}\right), \qquad \nu_k \sim \mathrm{IG}\!\left(\frac{1}{2},1\right), \qquad k=1,\ldots,d, \] and \[ \tau^2 \mid \xi \sim \mathrm{IG}\!\left(\frac{1}{2},\frac{1}{\xi}\right), \qquad \xi \sim \mathrm{IG}\!\left(\frac{1}{2},1\right), \] where $\mathrm{IG}(a,b)$ denotes the inverse-gamma distribution with density proportional to $ x^{-(a+1)}\exp\!\left(-\frac{b}{x}\right), x>0. $ Let $\boldsymbol{\nu}=(\nu_1,\ldots,\nu_d)^\top$. The joint posterior distribution of the latent states $\boldsymbol\theta_{1:T}$ and shrinkage parameters is given by
Algorithm (\oldref{alg:mwg_rw_hs_income_model}) in Section (\oldref{sec:additionaldetailshorseshoe}) of the online supplement summarises the MwG sampler for the random walk income model with horseshoe shrinkage priors. The two dynamic income models differ only in how the temporal innovation variances are modelled. In the RW model, each latent state has its own innovation variance $\sigma_k^2$. In the RW-HS specification, these variance terms are replaced by the product $\tau^2\lambda_k^2$, where $\tau^2$ controls the overall degree of smoothing and $\lambda_k^2$ allows each latent state component to depart from this global level.
This section summarises the inequality and poverty indices used in our empirical analysis. Given a parametric income model with density $p_{Y_t}(\cdot\mid\boldsymbol{\phi_t})$ and distribution function ${F_{Y_t}(\cdot\mid\boldsymbol{\phi_t})}$ at time $t$, we evaluate each index as a functional of $\boldsymbol{\phi}_t$ and a poverty line $z>0$. Let $\mu_t$ denotes the mean income, assumed finite, and $L_t(u;\boldsymbol{\phi}_t)$ denotes the Lorenz curve at time $t$, defined as the cumulative share of total income held by the poorest $u\in[0,1]$ proportion of the population.
The most widely used inequality measure is the Gini coefficient, $G_t$, defined as twice the area between the Lorenz curve and the line of equality. It ranges from $0$ to $1$, and is invariant to scale, in the sense that multiplying all incomes by a positive constant does not change $G_t$. The Gini index at time $t$ is given by
An equivalent expression, convenient for parametric modelling, is
see, for example, Gastwirth1971. In practice, (ref) is useful because closed-form expressions for $F_{Y_t}(\cdot\mid\boldsymbol{\phi})$ and for moment distribution functions often yield analytic or numerically stable evaluations of $G_t$ under parametric models.
We now briefly describe well-known poverty measures. Let $z>0$ denote a fixed poverty line, for example an absolute poverty threshold or a fraction of median income. We consider poverty measures that quantify complementary aspects of poverty: {incidence} (how many are poor), and {depth} (how far below the poverty line the poor are on average). These measures are standard in the poverty measurement literature; see, for example, FosterGreerThorbecke1984. The headcount ratio ($\mathrm{HC}_t$) is simply the proportion of the population with income below $z$ at time $t$:
The headcount ratio is easy to interpret but ignores the depth of poverty: it does not change when incomes of the poor fall further below $z$, provided they remain below the poverty line. The poverty gap ($\mathrm{PG}_t$) measures the average proportional shortfall from the poverty line at time $t$:
Unlike \(\mathrm{HC}_t\), the poverty gap accounts for how far incomes fall below the poverty line and can be interpreted as the minimum share of total income needed to raise all poor individuals to \(z\). More generally, the Foster--Greer--Thorbecke (FGT) class at time $t$ is defined by
Thus, $\mathrm{FGT}0_t=\mathrm{HC}_t$ and $\mathrm{FGT}1_t=\mathrm{PG}_t$. In this paper, we estimate the indices by Monte Carlo simulation: We draw a large sample from the fitted income distribution and compute the corresponding sample-based estimates.
This section discusses Lorenz dominance, generalised Lorenz dominance, and first-order stochastic dominance. Consider two income distributions at time $t$ indexed by parameter vectors $\boldsymbol{\phi}_{A,t}$ and $\boldsymbol{\phi}_{B,t}$, with corresponding cumulative distribution functions $F_{Y_{A,t}}(\cdot\mid\boldsymbol{\phi}_{A,t})$ and $F_{Y_{B,t}}(\cdot\mid\boldsymbol{\phi}_{B,t})$, means $\mu_{A,t}$ and $\mu_{B,t}$, and Lorenz curves $L_{A,t}(\cdot;\boldsymbol{\phi}_{A,t})$ and $L_{B,t}(\cdot;\boldsymbol{\phi}_{B,t})$. We state the dominance conditions using the population share $u\in[0,1]$. Distribution $A$ Lorenz (LD) dominates distribution $B$ at time $t$ if
Distribution $A$ generalised Lorenz (GLD) dominates distribution $B$ at time $t$ if
Distribution $A$ first-order (FSD) stochastically dominates distribution $B$ at time $t$ if
We follow lander2020bayesian and gunawan2021posterior to compute posterior probabilities of Lorenz and stochastic dominance. The posterior probabilities of Lorenz and stochastic dominance are estimated by counting the proportion of MCMC draws for which the estimated distributional functions satisfy the relevant dominance conditions. Since income is a continuous variable, when calculating the proportion of MCMC draws that satisfy the relevant dominance conditions for all $u\in [0,1]$, the best we can do is to check the conditions for a finite grid set of points. Accordingly, the resulting posterior probabilities should be interpreted relative to the adopted grid resolution and the total number of MCMC draws used in the analysis. In this paper, we see Section (\oldref{sec:additionaldetailLorenzstochastic}) for further details. In practice, this approach is analogous to widely used frequentist procedures, such as the tests of barrett2003consistent for first-order stochastic and generalized Lorenz dominance and barrett2014consistent for Lorenz dominance, which approximate continuum-based conditions by evaluating them over a discrete set of support points.
We also plot probability curves, which give the posterior probability of pointwise dominance at each population share $u$. Over any range of $u$, the posterior probability of dominance on that range can be no larger than the minimum value of the probability curve within the range. This makes the probability curve a useful device for identifying which parts of the distribution, such as the tails or the middle, drive the overall dominance probability. In particular, if dominance is largely determined by tail behaviour, one can assess sensitivity by omitting extreme values of $u$. Likewise, if interest focuses on a particular segment of the population, such as the poor, one can examine how the dominance probability changes when attention is restricted to the corresponding range of $u$.
This section presents a simulation study designed to assess how borrowing strength over time improves inference for income distributions modelled by the Dagum distribution Dagum1977. We generate annual incomes for $T=25$ years with $n=250$ observations per year, and we compare three models: an {independent} approach (ind) that fits a Dagum distribution separately for each year, and a {random-walk} (RW) approach in which transformed Dagum parameters evolve over time through a latent process following random walk model, and a random walk approach with horseshoe priors (RW-HS). All models are estimated in a fully Bayesian framework using the Markov chain Monte Carlo (MCMC) sampler discussed in Section (\oldref{sec:posteriorinference}).
Our evaluation assesses the ability of different models to estimate the true time-varying Dagum parameters, $(a_t,b_t,p_t)$, as well as a range of key distributional functionals. In particular, we examine posterior estimates of mean income $\mu_t$, the Gini coefficient $G_t$, and the $\mathrm{FGT}{0_t}$ and $\mathrm{FGT}{1_t}$ poverty indices, for $t=1,...,T$. We also compare the estimated densities, cumulative distribution functions, generalised Lorenz curves, and Lorenz curves over time obtained from posterior draws, with posterior credible bands reported for each curve.
We now discuss the data-generating process. We simulate data from the random-walk model specification. The log of the Dagum parameters follow a Gaussian random walk, producing gradual year-to-year changes in the income distribution. For $t=1,\ldots,T$ and $i=1,\ldots,n$, we assume $ y_{t,i}\mid (a_t,b_t,p_t) \overset{\text{iid}}{\sim} \mathrm{Dagum}(a_t,b_t,p_t). $ The initial values are set to $ a_1=3.54, b_1=329.58, p_1=0.61, $ and the corresponding latent log-parameters are defined as $ \tilde a_t=\log a_t, \tilde b_t=\log b_t, \tilde p_t=\log p_t. $ For $t=2,\ldots,T$, these latent log-parameters evolve according to $ \tilde a_t=\tilde a_{t-1}+\sigma_a \varepsilon_{a,t}, \tilde b_t=\tilde b_{t-1}+\sigma_b \varepsilon_{b,t}, \tilde p_t=\tilde p_{t-1}+\sigma_p \varepsilon_{p,t}, $ where $ \varepsilon_{a,t},\varepsilon_{b,t},\varepsilon_{p,t} \overset{\mathrm{iid}}{\sim}\mathcal N(0,1), $ and $ \sigma_a=\sigma_b=\sigma_p=0.02.$
Figure (\oldref{fig:Figure_param_dagum_sim1}) compares the log of the true parameters of the Dagum distribution with MCMC-based posterior summaries under the independent year-by-year model (ind), the random walk (RW) model, and the random walk model with horseshoe priors (RW-HS). For each parameter, the center curve represents the posterior mean, and the outer curves give the 95% credible intervals. The RW model tracks the true trajectories closely: it reproduces the gradual decline in $\log(b_t)$ and the mild fluctuations in $\log(a_t)$ and $\log(p_t)$, while delivering markedly tighter credible bands. In contrast, the independent model yields highly volatile year to year estimates and substantially wider uncertainty bands, most notably for $\log(p_t)$. The RW-HS model produces slightly smoother parameter trajectories than the RW model. Overall, the figure highlights the efficiency gains from modelling temporal dependence in the parameters.
Figure (\oldref{fig:Figure_dagum_mean_trajectory_sim1}) shows the true and estimated mean income, Gini coefficients, FGT0, and FGT1 indices over time obtained using the three models for the simulated dataset. The true mean income exhibits a clear downward trajectory over time, with only modest local fluctuations. The RW and RW-HS models track this gradual decline closely: their posterior means follow the true trajectory, and their pointwise credible bands are tight, with the RW-HS model producing slightly narrower intervals. In contrast, the independent model produces much noisier mean estimates and substantially wider credible bands, with several large year-to-year swings that are not present in the true trajectory. A similar pattern is observed for the headcount ratio (FGT0). The true FGT0 increases from the early years to the middle of the sample and then fluctuates only mildly thereafter. The RW and RW-HS models capture this smooth evolution with narrow credible intervals, whereas the independent model exhibits erratic movements and inflated uncertainty, leading to substantial overestimation or underestimation of the true values in some years.
For inequality, the true Gini coefficient trajectory remains relatively stable over the sample, with only small changes around its long-run level. The RW and RW-HS models again provide a close match to the true series, with smooth posterior means and narrow credible bands. By comparison, the independent model yields highly volatile Gini estimates, including pronounced spikes and dips (with wide bands).
The PDF and CDF plots in Figures (\oldref{fig:Figure_PDF_sim1}) and (\oldref{fig:Figure_CDF_sim1}) in Section (\oldref{sec:additionalfiguresimulationstudy}) of the online supplement (with 99% credible intervals) show that all models capture the shape of the true density and CDF. The independent model produces wide bands and large variations in peak height and location, indicating that year-by-year estimation with $n=250$ can translate sampling noise into spurious changes in the density and CDF. The RW and RW-HS models yield narrower credible bands and posterior means that are closer to the true density curve. The RW-HS model provides slightly tighter credible intervals than the RW model.
Figures (\oldref{fig:Figure_GLC_sim1_t25_099}) and (\oldref{fig:Figure_LC_sim1_t25_099}) in Section (\oldref{sec:additionalfiguresimulationstudy}) of the online supplement show the GLCs and LCs for time period $t=25$, obtained under the independent, RW, and RW-HS models. For the GLCs, uncertainty is smallest at low population shares and increases toward the upper tail; accordingly, the independent model exhibits a wider 99% credible band near the top of the distribution. The RW and RW-HS models deliver markedly tighter bands over the entire curve, especially for the upper deciles, implying more precise inference. The same conclusion holds for the LCs at $t=25$: the independent fit implies much greater uncertainty about income shares, particularly beyond the median and into the top quantiles, whereas the RW and RW-HS models produce narrower credible regions and smoother implied inequality profiles.
Tables (\oldref{tab:tableprobsim1}) and (\oldref{tab:tableprobsim1_poorest}) show that the estimated dominance probabilities can differ markedly across the independent model, the random walk model, and the random walk model with horseshoe shrinkage priors, with the two dynamic specifications generally producing much more similar results to each other than to the independent model. For the full population, the RW and RW-HS models tend to assign higher probabilities than the independent model for Lorenz dominance (LD), especially in the 2005--2001, 2015--2010, and 2025--2015 comparisons. For example, the LD probability for 2025 relative to 2015 increases from $0.2730$ under the independent model to $0.5667$ under RW and $0.5574$ under RW-HS, indicating substantially stronger evidence of inequality improvement once temporal smoothing is imposed. A similar pattern appears for 2005 relative to 2001 and 2015 relative to 2010, where the dynamic models roughly double the LD probabilities relative to the independent specification. By contrast, for first-order stochastic dominance (FSD) and generalised Lorenz dominance (GLD), the most striking discrepancy occurs for the 2010--2005 comparison. Here the independent model suggests moderate evidence of welfare improvement, with probabilities $0.1211$ for FSD and $0.1974$ for GLD, whereas the RW and RW-HS models reduce these probabilities to values close to zero. This indicates that the independent model may overstate welfare improvement. For the remaining comparisons, FSD and GLD probabilities are generally low under all three models, although the RW and RW-HS specifications often produce slightly larger probabilities than the independent model for 2005--2001 and 2025--2015.
The differences across models are even more pronounced when attention is restricted to the poorest 10% of the population. In this case, the independent model often yields substantially larger probabilities of FSD and GLD than the dynamic models, particularly for the 2010--2005 and 2025--2015 comparisons. For example, for 2010 relative to 2005, the independent model reports probabilities of $0.4511$ for FSD and $0.6004$ for GLD, whereas the corresponding probabilities under RW and RW-HS fall dramatically to around $0.02$ or lower. Likewise, for 2025 relative to 2015, the independent model gives much stronger evidence of dominance at the lower end of the distribution than either dynamic specification. In contrast, the RW and RW-HS give substantially smaller dominance probabilities. However, the dynamic models give larger probabilities of Lorenz dominance relative to the independent model in several cases, such as 2005--2001 and 2015--2010. The RW and RW-HS results are very close throughout.
Figures (\oldref{fig:yFSDx_prob_sim1})--(\oldref{fig:yLDx_prob_sim1}) show that the probability curves differ systematically across the independent, RW, and RW-HS specifications, and these differences help explain the dominance probabilities reported in Tables (\oldref{tab:tableprobsim1}) and (\oldref{tab:tableprobsim1_poorest}). For FSD and GLD, the independent model generally produces much more irregular and more extreme pointwise probability curves, with pronounced U-shaped or hump-shaped patterns across the support. This is especially evident for the 2010--2005 and 2025--2015 comparisons, where the independent model assigns relatively high pointwise dominance probabilities near the tails but much lower probabilities in the middle of the distribution. By contrast, the RW and RW-HS curves are markedly smoother and less erratic, reflecting the stabilising effect of temporal pooling. In several cases, such as 2010--2005, the dynamic models produce probability curves that remain close to zero over most of the support, indicating that once year-to-year noise is smoothed out there is little evidence of dominance. More generally, these figures show that the independent model is much more sensitive to local fluctuations in the estimated income distributions, whereas the RW and RW-HS models deliver more coherent probability profiles over the whole range of population shares.
The LD curves reveal a somewhat different pattern. Here, the RW and RW-HS models often lie above the independent model over large parts of the support, particularly for the 2005--2001, 2015--2010, and 2025--2015 comparisons, which is consistent with their substantially higher posterior probabilities of LD in Table (\oldref{tab:tableprobsim1}). The RW and RW-HS curves are again very similar, showing that the horseshoe prior mainly refines the degree of smoothing rather than changing the substantive conclusions. At the same time, the figures also illustrate why quite high pointwise probabilities do not necessarily translate into high overall dominance probabilities: dominance must hold simultaneously for all evaluation points, so even a relatively small dip in the curve can substantially reduce the joint posterior probability.
In summary, the simulation results show that the RW and RW-HS models substantially outperform the independent model. Both dynamic specifications track the true parameter and distributional trajectories much more closely, with smoother estimates and much narrower credible intervals, while the independent model produces noisy estimates and inflated uncertainty. These improvements also carry over to the dominance analysis, where the RW and RW-HS models yield more stable and coherent probability profiles. Overall, modelling temporal dependence yields more accurate and reliable inference, with RW-HS providing an additional gain over RW.
Section (\oldref{sec:RealDataApplication}) briefly describes HILDA data. Section (\oldref{sec:aboriginalpopulationsubgroup}) discusses the empirical results for the Aboriginal population subgroup.
We use HILDA data from 2001 to 2021 to study the income distributions over time of Aboriginal subpopulations in Section (\oldref{sec:aboriginalpopulationsubgroup}) and residents of the Australian Capital Territory (ACT) in Section (\oldref{sec:ACTpopulationsubgroup}) of the online supplement. The welfare measure considered is household disposable income converted to a per-individual basis. Net disposable household income is obtained by subtracting “total disposable income negative per household” from “total disposable income positive per household”. Following SilaDugain2019, equivalised income is constructed by dividing net disposable household income by the square root of household size, and this equivalised value is then assigned to each household member. The analysis is restricted to individuals aged 15 years. However, children aged below 15 are still included in the household size used to calculate equivalised income. Values were adjusted using the Consumer Price Index, treating 2021/2022 as the base.
This section presents the empirical application of our modelling framework to Aboriginal population subgroups. Using HILDA data from 2001 to 2021, we compare a range of parametric income models to determine which specification best captures the evolution of the income distributions over time. Specifically, we consider the Dagum, Singh--Maddala, Beta 2, and GB2 distributions, described in Section (\oldref{sec:incomedistributions}) of the online supplement, under an independent model, a random walk model, and a random walk model with horseshoe shrinkage priors.
Table (\oldref{tab:predictivescoreabo}) reports the log predictive scores GneitingRaftery2007 obtained from 10-fold cross-validation for the Aboriginal population subgroup, where larger values indicate better out-of-sample predictive performance. For each year, the observations are randomly divided into ten approximately equal-sized folds. For a given fold, the model is fitted using the observations in the remaining nine folds from each year, and the fitted model is then used to evaluate the predictive density of the held-out observations. This procedure is repeated ten times so that each fold is used once as the validation sample. The log predictive score is then computed by averaging the log of the predictive densities of all held-out observations across all years. Several conclusions emerge from the table. First, allowing the income distribution parameters to evolve over time generally improves predictive accuracy relative to the independent model. This is evident for the Dagum, Beta 2, and GB2 specifications, for which both the RW and RW-HS models achieve higher log predictive scores than the corresponding independent model.
Second, among the four candidate income distributions, the GB2 model performs best overall, with the RW-HS specification attaining the highest log predictive score of $-4009.22$, followed very closely by the RW version at $-4009.32$. This suggests that the random walk GB2 model provides the most flexible and accurate representation of the income distributions for this subgroup, while the addition of horseshoe shrinkage yields a further small improvement. Third, the Dagum model also performs competitively, with RW and RW-HS log predictive scores of $-4009.70$ and $-4009.80$, respectively, although both are slightly worse than the corresponding GB2 specifications.
In contrast, the Beta 2 model performs substantially worse than the other candidates under all three model structures, with log predictive scores around $-4032$ to $-4035$, indicating poor predictive fit. Finally, the Singh--Maddala model shows almost no difference across the independent, RW, and RW-HS specifications, suggesting that introducing temporal dependence or shrinkage provides little benefit for this distribution in the present application. Overall, the cross-validation results support selecting the RW-HS GB2 model as the preferred specification for the Aboriginal population subgroup, and it is therefore used as the main model in the discussion that follows.
Figures (\oldref{fig:Figure_param_GB2_abo}) and (\oldref{fig:Figure_GB2_mean_gini_FGT0_FGT1_abo}) present the estimated GB2 parameters and the corresponding welfare measures over time for the Aboriginal population subgroup under the independent model (ind) and the random walk model with horseshoe shrinkage priors (RW-HS). The most notable difference between the two specifications is that the RW-HS model yields considerably smoother and more stable temporal trajectories, whereas the independent model exhibits substantial year-to-year fluctuations, particularly in the estimates of \(\log(a)\) and \(\log(p)\). Despite these differences in smoothness, both models imply similar broad distributional trends. In particular, mean income increases over time, especially in the later years of the sample, while the Gini coefficient indicates a moderate rise in inequality in the earlier period, followed by relative stability thereafter. The FGT0 and FGT1 measures both display a general downward trend, pointing to declines in the incidence and intensity of poverty over the study period for the Aboriginal population subgroup.
Figures (\oldref{fig:Figure_PDF_GB2_abo})--(\oldref{fig:Figure_LC_GB2_abo}) in Section (\oldref{sec:additionalfiguresaboriginalsubgroup}) of the online supplement present the posterior means and 99% credible intervals for the fitted GB2 income PDFs, CDFs, GLCs, and LCs for the Aboriginal population subgroup at time \(t=21\) under the independent model and the random walk model with horseshoe priors. The posterior mean CDFs, LCs, GLCs, and PDFs under the independent model and the random walk model with horseshoe priors are very close to one another, indicating that both models deliver essentially the same overall distributional picture. The differences in posterior means are generally small rather than substantial: the RW-HS LC lies slightly above that of the independent model over most population shares, suggesting marginally lower inequality, while the corresponding GLC is slightly lower, reflecting a somewhat smaller fitted mean income. Likewise, the posterior mean CDFs and PDFs are very similar across the support, with only minor deviations in the lower and middle parts of the income distribution. Importantly, the RW-HS model tends to produce narrower 99% credible intervals than the independent model, indicating greater estimation precision and a more stable characterisation of uncertainty.
Tables (\oldref{tab:FSDGLDLDOVERALL}) and (\oldref{tab:FSDGLDLDPOOREST}) show that the posterior probabilities of Lorenz and stochastic dominance can differ substantially between the independent model and the random walk model with horseshoe priors, especially for intermediate year-to-year comparisons and for the poorest 10% of the population. For FSD, both models agree that there is essentially no evidence that the 2005 income distribution dominates that of 2001, either overall or for the poorest 10%, and both models also agree that the strongest welfare improvement occurs from 2015 to 2021, with posterior probabilities close to one overall and effectively one for the poorest 10%. The main differences arise in the middle of the sample. From 2005 to 2010, the independent model suggests moderate evidence of FSD overall and very strong evidence for the poorest 10%, whereas RW-HS yields much weaker support. By contrast, from 2010 to 2015, the independent model assigns essentially zero probability to FSD, while RW-HS gives noticeably larger probabilities, particularly for the poorest 10%, indicating some improvement at the lower end of the distribution.
The GLD results tell a similar story: there is almost no evidence of welfare improvement from 2001 to 2005, very strong evidence of improvement from 2015 to 2021, and markedly different conclusions across models for 2005--2010 and 2010--2015, with the independent model favouring improvement in the former period and RW-HS giving relatively more support in the latter. For LD, which focuses on inequality comparisons, the evidence is generally weaker at the overall population level, suggesting that changes in welfare were not driven primarily by uniform reductions in inequality across the whole distribution. However, the poorest 10% display a more nuanced pattern: the independent model suggests strong inequality improvement from 2005 to 2010, RW-HS instead gives stronger support for inequality improvement from 2010 to 2015, and both models indicate very strong inequality improvement for the poorest 10% from 2015 to 2021. Overall, these results suggest that the most robust welfare gain for the Aboriginal population subgroup occurs between 2015 and 2021.
The probability curves in Figures (\oldref{fig:yFSDx_prob_abo}) to (\oldref{fig:yLDx_prob_abo}) in Section (\oldref{sec:additionalfiguresaboriginalsubgroup}) of the online supplement reinforce these differences between the independent model and RW-HS by showing that the two specifications can imply very different pointwise dominance behaviour, even when the resulting overall dominance probabilities are similar in broad direction. Across the FSD, GLD, and LD panels, the RW-HS curves are generally smoother and less extreme, whereas the independent model often produces sharply varying, and hump-shaped profiles. This is particularly evident for the 2015--2010 comparison, where the independent model yields highly uneven probability curves across the support, while RW-HS gives flatter and more regular profiles. For 2021 versus 2015, by contrast, both models produce probability curves that are close to one throughout most of the support for FSD and GLD, consistent with the strong dominance probabilities reported in the tables. Another important feature is that pointwise dominance probabilities can be high over large parts of the support while the joint dominance probability remains low, because dominance must hold simultaneously at all evaluation points. This is especially clear for the 2005--2001 and several LD comparisons, where the curves may be moderately large over part of the domain but still fail to imply a large overall probability of dominance.
Figures (\oldref{fig:pred_density_Abo})--(\oldref{fig:pred_LC_Abo}) in Section (\oldref{sec:additionalfiguresaboriginalsubgroup}) of the online supplement present the posterior predictive distributions for the Aboriginal population subgroup in 2022 and 2025, obtained by projecting the RW-HS GB2 model beyond the observed 2001--2021 period. Figure (\oldref{fig:pred_density_Abo}) shows that the 2025 predictive density is more dispersed than that for 2022, with a noticeably wider range of plausible incomes and a more pronounced upper tail, indicating greater uncertainty about the future shape of the income distribution. This pattern is also reflected in Figure (\oldref{fig:pred_CDF_Abo}), which shows the predictive CDFs for 2022 and 2025. Figures (\oldref{fig:pred_GLC_Abo}) and (\oldref{fig:pred_LC_Abo}) suggest that the predicted generalised Lorenz and Lorenz curves for 2002 closely resemble those for 2025. Importantly, however, the prediction intervals are clearly wider in 2025 than in 2022 across all four curves, showing that forecast uncertainty accumulates substantially as the prediction horizon moves further away from the sample used for estimation.
Figures (\oldref{fig:pred_means_Abo})--(\oldref{fig:pred_FGT1_Abo}) in Section (\oldref{sec:additionalfiguresaboriginalsubgroup}) of the online supplement summarise the predictive distributions of key welfare measures for the Aboriginal population subgroup over the out-of-sample period 2022--2025 under the RW-HS GB2 model fitted to the 2001--2021 data. Figure (\oldref{fig:pred_means_Abo}) shows that the predicted mean income remains centred in a broadly similar range across the forecast horizon, although the predictive density becomes progressively flatter and more dispersed from 2022 to 2025, indicating that uncertainty about future mean income increases substantially for the later years. A similar pattern is evident in Figure (\oldref{fig:pred_GINI_Abo}), where the predictive densities for the Gini coefficient remain concentrated around comparable values, suggesting no dramatic change in relative inequality, but the later-year densities are clearly wider and exhibit more tail mass, so any apparent movement should be interpreted cautiously. The poverty measures in Figures (\oldref{fig:pred_FGT0_Abo}) and (\oldref{fig:pred_FGT1_Abo}) suggest a modest tendency towards lower poverty by 2025, as the predictive mass shifts slightly toward smaller values for both the headcount ratio and the poverty gap. However, these improvements are accompanied by much greater dispersion and longer right tails in the later years, especially for 2024 and 2025, reflecting the build-up of forecast uncertainty as the prediction horizon extends further beyond the observed sample. Overall, the model points to broadly stable or slightly improving welfare outcomes for the Aboriginal subgroup, but the substantially wider predictive densities in the later years make clear that these longer-horizon forecasts are much less precise than those for 2022.
This paper develops a flexible Bayesian time-varying parametric model for income distributions in which the parameters of income distributions are allowed to evolve over time, rather than being estimated separately and independently for each year. By embedding the income distribution parameters in a latent random-walk state process, the proposed framework borrows strength across adjacent years and thereby produces more stable and coherent inference for inequality, poverty, and welfare comparisons. This is especially important for population subgroups with relatively small numbers of observations, where independent year-by-year estimation can generate noisy parameter paths, wider credible intervals, and spurious year-to-year variation in derived welfare measures. In contrast, the proposed model preserves the interpretability and tractability of parametric income modelling while substantially improving estimation precision and allowing coherent prediction of future income distributions and related welfare, inequality, and poverty summaries. The simulation results reinforce these advantages, showing that the dynamic specification tracks the underlying time-varying parameters more closely and delivers tighter uncertainty quantification than the corresponding independent income models, particularly when the true distribution evolves smoothly over time. A further contribution of the paper is the use of horseshoe shrinkage priors within the dynamic income modelling framework. The horseshoe prior provides smoother temporal trajectories for the income distribution parameters and for the associated inequality and poverty measures.
In the empirical application, this benefit is clearly visible for both the Aboriginal population subgroup and the Australian Capital Territory (ACT) subgroup, where the random walk GB2 model with horseshoe priors is selected by 10-fold cross-validation as the preferred specification. For the Aboriginal subgroup, the results indicate rising mean income over time, a moderate increase in inequality in the earlier years followed by relative stability, and an overall decline in the incidence and intensity of poverty. The dominance analysis further shows that the strongest and most robust welfare improvement occurs between 2015 and 2021, particularly for the poorest 10% of the population. For the ACT subgroup, the preferred model similarly yields smoother parameter estimates and more stable welfare summaries, with mean income generally increasing over time, inequality rising in the earlier part of the sample and then stabilising or declining slightly, and poverty measures exhibiting an overall downward trend. The dominance results suggest that the clearest welfare gains for the ACT subgroup also occur in the later part of the sample, especially at the lower end of the income distribution. Overall, the findings demonstrate that the proposed time-varying income model with horseshoe shrinkage priors provides a practical and useful alternative to independent income models for analysing changes in income distributions over time. For longer time series, this framework could be extended further by using dynamic shrinkage priors, which can capture both smooth parameter evolution and occasional abrupt changes or jumps in the income distribution parameters kowal2019dynamic,knaus2023dynamic.
\setcounter{page}{1} \setcounter{section}{0} \setcounter{equation}{0} \setcounter{algorithm}{0} \setcounter{table}{0} \setcounter{figure}{0}