---
title: Modeling Income and Wealth Distributions I
author: 36-313, Fall 2022
date: 6 September 2022 (Lecture 3)
output:
  pdf_document:
    toc: true
    number_sections: true
---

\newcommand{\Prob}[1]{\mathbb{P}\left( #1 \right) }
\newcommand{\Expect}[1]{\mathbb{E}\left( #1 \right) }
\newcommand{\Var}[1]{\mathrm{Var}\left( #1 \right) }
\newcommand{\myexp}[1]{\exp{\left\{ #1 \right\}}}
\newcommand{\xmin}{x_{\min}}

```{r, include=FALSE}
# Set up global options for displaying code and so forth
library(knitr)
opts_chunk$set(size="small", 
               tidy=TRUE, tidy.opts=list(comment=FALSE),
               cache=TRUE, cache.comments=FALSE, autodep=TRUE,
               warning=FALSE, message=FALSE, error=FALSE, echo=FALSE)
```

In this lecture, we're going to look at _models_ for the distribution of income
and wealth.  These aren't facts about they world; they're ideas we've invented
to help organize and make sense of those facts.  Models are simplified,
stylized pictures which deliberately suppress some details in order to
highlight what are (hopefully) larger and more important patterns.  To decide
which stylized picture to draw, though, we need to pay some attention to
the facts.

\clearpage

# Reminders of what income distributions look like

```{r, include=FALSE}
# Repeated from lecture 2

# Data download and re-shaping

# The WID data is organized by percentile _ranges_, so the names we want always
# need to give lower and upper limits of the range.  But we just care about
# the threshold which puts you at a given percentile, so we can always take
# the upper limit to be 100.
percentile.levels <- c(1:99, seq(from=99.1, to=99.9, by=0.1),
            seq(from=99.91, to=99.99, by=0.01),
            seq(from=99.991, to=99.999, by=0.001))
percentile.names <- paste("p", percentile.levels, "p100", sep="")

# Download thresholds (t) for pretax  income (ptinc) in the US and selected
# other countries over all
# available years, at every available percentile, with a special focus on the
# very top.
  # Limiting to adults (but not to special age ranges of adults)
  # Assigning equal shares to all members of households
library(wid)
some.countries <- c("US","DE","FI","FR","GB","JP", "BR","CN","ZA", "SU")
ctryinc <- download_wid(indicators="tptinc", areas=some.countries,
                         years=1948:2019,
                         ages=992, # Adults only (but all adults)
                         perc=percentile.names,
                         pop="j")

# Note: at this point, within each year the variables are apparently arranged
# in _alphabetical_ order, so that for instance "p99.9p100" comes _before_
# the numerically-lower "p9p100".  I'll fix this later.


# WID returns the data in a "long" format, one row per year-country-variable
# combination
# It'll be more convenient to reshape this to a "wide" format, with one
# row per year-country pair, and multiple columns for multiple values
# The following, taken from [https://stackoverflow.com/questions/18047687/reshape-long-data-frame-into-a-wide-structure-using-plyr], does this, using
# the setNames() function (which was new to me)
library(plyr)
ctryinc.wide <- ddply(ctryinc, .(year, country),
                       .fun=function(t) {
                           setNames(object=t$value, nm=t$percentile)
                       })


# There is now one column per variable, but in the messed-up alphanetized
# order, so let's put things right by re-arranging the columns
ctryinc.wide <- ctryinc.wide[,c("year", "country", percentile.names)]

# Now fix the column names to be less horrible
  # I should really use search-and-replace with pattern matching here,
  # but it's too much work to remember how that goes in R
colnames(ctryinc.wide)[-(1:2)] <- paste("P", percentile.levels, sep="_")

# The reshaping looks good, so change the name
ctryinc <- ctryinc.wide
```

Here, [as a
reminder](http://www.stat.cmu.edu/~cshalizi/ineq/21/lectures/01/lecture-01.pdf),
is an approximation to the distribution of individual income for the US in
2019.

```{r, zoomedinpdf}
usinc <- as.numeric(ctryinc[ctryinc$year==2019 & ctryinc$country=="US", -(1:2)])
plot(x=tail(usinc, -1),
     y=diff(percentile.levels/100)/diff(usinc),
     xlim=c(-2000, 4e5), type="l",
     xlab="Taxable income ($)", ylab="Probability density",
     main="Taxable income per individual, 2019")
abline(h=0, col="grey")
axis(side=3, at=tail(usinc,-1), labels=tail(percentile.levels,-1),
     cex.axis=0.5)
```

If we zoom out to the full range, this looks pretty much like a spike,
because the range is so great:

```{r}
plot(x=tail(usinc, -1),
     y=diff(percentile.levels/100)/diff(usinc),
     type="l",
     xlab="Taxable income ($)", ylab="Probability density",
     main="Taxable income per individual, 2019")
abline(h=0, col="grey")
axis(side=3, at=tail(usinc,-1), labels=tail(percentile.levels,-1),
     cex.axis=0.5)
```

On the other hand, the little numbers on the upper edge of the plot show percentiles of income, so it's not that there's _no_ probability density out around an income of (say) 20 million dollars a year.  We can get a bit better sense of what's going on by using a logarithmic scale on both axes:

```{r, loglogpdf}
plot(x=tail(usinc, -1),
     y=diff(percentile.levels/100)/diff(usinc),
     xlim=c(100, max(usinc)), type="l", log="xy",
     xlab="Taxable income ($)", ylab="Probability density",
     main="Taxable income per individual, 2019")
abline(h=0, col="grey")
axis(side=3, at=tail(usinc,-1), labels=tail(percentile.levels,-1),
     cex.axis=0.5)
```

What should we take away from these plots?

1. There's fairly substantial probability density for a wide range around the
median income, maybe from an order of magnitude below it to an order of magnitude above it.  (Look at the close bunching of the quantile marks in the log-scale plot.)
2. The distribution is highly right-skewed: the mean is much larger than the median.
3. The distribution is heavy tailed: high quantiles are many orders of magnitude larger than the median, and very high quantiles are orders of magnitude larger than merely high quantiles.
4. At very high incomes, the log-log plot of density looks like a straight line.  This is what we'd see if the pdf were shrinking like some power of income, $f(x) \propto x^{-\alpha}$, since then $\log{f(x)} = -\alpha\log{x} + \text{(intercept)}$.

All of these are pretty common features of income and wealth distributions.
(See Complementary Problem 1.)  When we introduce models, it'd be nice if they
could account for, or at least reproduce, some of these features.  While
there isn't, unfortunately, one simple model which does all of this, there are
two simple models each of which captures _part_ of what's going on.

\clearpage

# The Log-Normal Distribution

You remember the Gaussian or "normal" distribution; it's conveniently
defined by its pdf,
\[
f(x) = \frac{1}{\sqrt{2\pi\sigma^2}}\myexp{-\frac{(x-\mu)^2}{2\sigma^2}}
\]
There are two adjustable knobs in this, the **parameters** $\mu$ and
$\sigma^2$.  We write this $X \sim \mathcal{N}(\mu, \sigma^2)$, where the
symbol $\sim$ is read "is distributed as".

Some calculus, which I won't reproduce, shows that with this pdf, $\Expect{X} =
\mu$ while $\Var{X} = \sigma^2$.  The highest probability
density, the mode, occurs when $x=\mu$.  As $x$ moves away from $\mu$ in either direction, the pdf falls off, because $(x-\mu)^2 \geq 0$.  Initially, the fall-off is slow, because if $|x-\mu|/\sigma < 1$, then $(x-\mu)^2/\sigma^2 < |x-\mu|/\sigma$.  (Incidentally, this algebra explains why "close to $\mu$" and "far from $\mu$" are measured with respect to the standard deviation $\sigma$.)  But conversely, once $|x-\mu|/\sigma$ becomes large, the pdf goes to zero very rapidly.  The familiar conclusion is that with high probability, a random $X$ drawn from this distribution will be found near $\mu$, meaning within $\pm$ a few multiples of $\sigma$ of $\mu$.

The log-normal distribution is just what it sounds like: $X$ has the
log-normal[^hyphen] distribution when $\log{X}$ has a normal or Gaussian
distribution[^logbase].  We write this as $X \sim \mathcal{LN}(\mu, \sigma^2)$,
with the understanding that now $\Expect{\log{X}} = \mu$, $\Var{\log{X}} =
\sigma^2$.  What does this mean in practice?

[^hyphen]: Some people write "log-normal", others write "lognormal"; no difference.  I tend to use the hyphen because it makes my spell-checker happier.

[^logbase]: It really doesn't matter what "base" of logarithm we're using, whether $\log$ is really $\log_{10}$ or $\log_2$ or $\ln$.  But some of the calculus is easier if we use natural logs, so if you want to read $\ln$ for $\log$ throughout, that's fine.

## Some properties of log-normals

First, $\log$ is a strictly increasing function, so
\[
\Prob{X \leq b} = \Prob{\log{X} \leq \log{b}}
\]
and thus
\[
\Prob{a < X \leq b} = \Prob{\log{a} < \log{X} \log{b}}
\]
One immediate implication is that quantiles of $X$ are just the exponentials of
the quantiles of $\log{X}$.  In particular, the median of $X$ is $e^{\mu}$.
(Why does that follow?)

Another implication of this basic relationship is the following.  We know that it's very improbable for $\log{X}$ to be more than a few multiples of $\sigma$ away from $\mu$.  Let's write this more formally:
\[
\Prob{\mu - k \sigma < \log{X} \leq \mu + k \sigma} \approx 1
\]
where $k$ is just a convenient algebraic symbol for "a few".  But this implies
\[
\Prob{e^{\mu} (e^{-\sigma})^k < X \leq e^{\mu} (e^{\sigma})^k} \approx 1
\]
So it's very likely that $X$ is within a _factor_ of a few powers of $e^{\sigma}$ of $e^{\mu}$.  The range probable range of $X$, in other words, is multiplicative, not additive.  (This is already promising in view of the empirical pdf of income.)

We can quickly find the exact pdf of the log-normal
distribution.  Remember the definitions: the CDF is $F(a) \equiv \Prob{X \leq a}$, and the pdf is the derivative of the CDF, $f(x) = \left.\frac{dF}{da}\right|_{a=x}$.  So when $X$ is log-normal,
\begin{eqnarray}
f(x) & = & \left. \frac{d}{da}\Prob{X \leq a} \right|_{a=x}\\
& = & \left. \frac{d}{da}\Prob{\log{X} \leq \log{a}} \right|_{a=x}\\
& = & \frac{1}{\sqrt{2\pi\sigma^2}}\myexp{-\frac{(\log{x}-\mu)^2}{2\sigma^2}} \frac{d\log{x}}{dx}\\
& = & \frac{1}{x}\frac{1}{\sqrt{2\pi\sigma^2}}\myexp{-\frac{(\log{x}-\mu)^2}{2\sigma^2}}
\end{eqnarray}

Now that we have the pdf, we can plot it.  Here's the pdf of the "standard"
log-normal, $\mathcal{LN}(0, 1)$:

```{r}
curve(dlnorm(x), from=0, to=10, xlab="x", ylab="pdf")
```

Notice that while the median is $e^0=1$, the mode is clearly lower than that.
Since we know the pdf, we can find the mode exactly by solving $\frac{df}{dx} =
0$ (Complementary Problem 2a), with the result that the mode is exactly
$e^{\mu - \sigma^2}$.

On the other hand, mean of the log-normal is higher than the median, $e^{\mu +
\sigma^2/2}$ (Complementary Problem 2b).  This reassures us that, as we'd
expect from the plot of the pdf, the log-normal distribution is right-skewed,
since $\sigma^2 > 0$.

## Why we care: log-normal versus income data

Here's our first plot of the pdf of income (as calculated from the income percentiles) again, with a log-normal distribution added in blue.

```{r}
usinc.median <- with(ctryinc,
                      ctryinc[country=="US" & year==2019, "P_50"])
# Get the mean US individual income from the same source
usinc.mean <- download_wid(indicators="aptinc", areas="US",
                         years=2019, perc="p0p100", pop="i")$value
```

```{r, eval=TRUE}
<<zoomedinpdf>>
curve(dlnorm(x, meanlog=log(usinc.median),
             sdlog=sqrt(2*log(usinc.mean/usinc.median))), col="blue",
      add=TRUE)
```

(This log-normal is fitted in exactly the way you're learning to do in HW 1.
We'll cover fitting methods in more detail next lecture.)  What you can see is
that while the log-normal distribution isn't quite right at very low incomes,
and in particular messes up how many people have incomes close to 0, from about
the 10th percentile of the income distribution onwards it's a close match
to the data, at least visually, through the 98th percentile and perhaps above.

Let's look at this same comparison in a log-log plot, which will emphasize
the right tail.

```{r}
<<loglogpdf>>
curve(dlnorm(x, meanlog=log(usinc.median),
             sdlog=sqrt(2*log(usinc.mean/usinc.median))), col="blue",
      add=TRUE)
```

Here we can see that there's a visually good match over most of the range of
incomes, from about the 10th percentile (or maybe a little lower?) to above the
99th percentile.  The fact that the log-normal distribution systematically
under-predicts how many people have (or report) very low incomes shows up
again.  What's new is being able to see what's happening in the right or upper
tail: the log-normal also under-predicts how many people have extremely high
incomes.  (Said another way, it predicts that the highest percentiles of income
should be much lower than they actually are.)  So while the log-normal is a
pretty good approximation to the distribution of income for most of the
population, "the rich are different than you and me".

\clearpage

# The Power-Law or Pareto Distribution

A good (but not perfect) model of the upper tail of income and wealth
distributions has actually come down to us from the first social scientist to
seriously investigate the issue, the pioneering economist and sociologist
Vilfredo Pareto, who first published his findings in the 1890s.  He used tax
records to examine the number of people with very high incomes or net worths in
a number of countries and cities for which he could find the data.  What he did
visually was plot the number of tax-payers whose income met or exceeded any
given level (e.g., $1000, $2000, ...) against that threshold level.  In other
words, he plotted $\Prob{X \geq x}$ against $x$.  In modern terminology, we'd
call this the **complementary cumulative distribution function** (CCDF) or
**upper CDF** or **survival function**, often written $\overline{F}(x)$.  When
Pareto made these plots with log scales on both the vertical and horizontal
axes, he found that the curve became a straight line, at least above some
minimum income or wealth.  This implied that $\overline{F}(x)$ was proportional
to some power of $x$, at least above that minimum level.

What we now call the **continuous power law distribution** or **Pareto distribution** has the following pdf:
\[
f(x) = \left\{ \begin{array}{cc} 0 & x < \xmin\\
\frac{\alpha-1}{\xmin} \left(\frac{x}{\xmin}\right)^{-\alpha} & x \geq \xmin
\end{array} \right.
\]
In words, the pdf is a power of $x$, hence the name "power law".

If you integrate this, you will find that, for $x \geq \xmin$,
\[
\Prob{X \geq x} = \left( \frac{x}{\xmin} \right)^{-\alpha+1}
\]
so the CCDF is _also_ a power of $x$.

At this point I need to say that there is a small but annoying discrepancy in
the literature.  Basically everyone writes $\alpha$ for the Pareto exponent,
but for some people, that's the exponent in the pdf, and for other people,
that's the exponent in the CCDF.  The two exponents always differ by exactly
one, so this isn't a big deal, but if you just want to borrow a formula, or
compare empirical results, it's annoying to have to constantly check which
convention is being used by a particular author.

## Some properties of the Pareto distribution

### Mode

The mode of the Pareto distribution is $\xmin$: the highest probability density is _always_ at the smallest possible value.

### Median and other quantiles

If we assume that there's no probability mass below $\xmin$, then it's easy
to find upper quantiles:

\begin{eqnarray}
\Prob{X \geq x_q} & = & q\\
\left(\frac{x_q}{\xmin}\right)^{-\alpha+1} & = & q\\
x_q & = & \xmin q^{1/(1-\alpha)}
\end{eqnarray}

Thus the ordinary, lower quantiles are also easy:
\begin{eqnarray}
\Prob{X \leq x_p} & = & p\\
x_p & = & \xmin (1-p)^{1/(1-\alpha)}
\end{eqnarray}

Matters are a little more complicated if we're just using the Pareto as a model
for the tail above $\xmin$.  _If_ $x_q \geq \xmin$, then
\begin{eqnarray}
\Prob{X \geq x_q} & = & \Prob{X \geq x_q|X \geq \xmin}\Prob{X \geq \xmin}\\
& = & \left(\frac{x_q}{\xmin}\right)^{-\alpha+1}\Prob{X \geq \xmin}
\end{eqnarray}
Setting this equal to $q$ and solving for $x_q$ gives
\begin{eqnarray}
x_q & = & \xmin (q/\Prob{X \geq \xmin})^{1/(1-\alpha)}
\end{eqnarray}

### Expectation and higher moments

The expected value of $X$ is easily calculated:
\begin{eqnarray}
\Expect{X| X \geq \xmin} & = & \int_{\xmin}^{\infty}{x \frac{\alpha-1}{\xmin} \left(\frac{x}{\xmin}\right)^{-\alpha} dx}\\
& = & \frac{\alpha-1}{\xmin} \xmin^{\alpha} \int_{\xmin}^{\infty}{x^{-\alpha+1}  dx}\\
& = & (\alpha-1)\xmin^{\alpha-1} \frac{1}{\alpha-2}\xmin^{-\alpha+2}\\
& = & \xmin \frac{\alpha-1}{\alpha-2}
\end{eqnarray}

Notice that as $\alpha \rightarrow 2$, this goes to infinity.  The practical
meaning of this is that when $\alpha \leq 2$, as we get larger and larger
samples, the sample mean will keep getting larger and larger without any limit.
This is _not_ how most random variables work, but it is the way of very
heavy-tailed distributions.

Higher moments are even less well-behaved:
\begin{eqnarray}
\Expect{X^k|X \geq \xmin} & = & \int_{\xmin}^{\infty}{x^k \frac{\alpha-1}{\xmin} \left(\frac{x}{\xmin}\right)^{-\alpha} dx}\\
& = & (\alpha-1)\xmin^{\alpha-1}\int_{\xmin}^{\infty}{x^{k-\alpha} dx}\\
& = & \frac{\alpha-1}{\alpha-k-1}\xmin^{\alpha-1}\xmin^{k+1-\alpha}\\
& = & \xmin^k \frac{\alpha-1}{\alpha-k-1}
\end{eqnarray}

Since the variance is $\Var{X} = \Expect{X^2} - (\Expect{X})^2$, it follows
that the variance is infinite if $\alpha \leq 3$.  Again, the practical meaning
is that if you keep drawing larger and larger samples from such a distribution,
the sample variance will grow without any limit.

## Why we care: Pareto distributions versus tail data

You will explore this in detail in Homework 2.  I will just say here that while
Pareto distributions aren't a _perfect_ model fit to the upper tail, they are
pretty good approximations.  You can, for instance, see this in the log-log
plots of the pdf I gave above --- the right tail is very close to a straight
line.

\clearpage

# Calculating measures of inequality from theoretical distributions

One use of these models is to simplify or "regularize" the calculation of
Lorenz curves, Gini indices, etc.

## Lorenz curves and Gini indices from the Pareto distribution

Let's do an example of calculating Lorenz curves and Gini coefficients
from a Pareto distribution.  To simplify matters, I will assume that there
is _no_ probability mass below $\xmin$, but you can re-do the calculations
without this assumption (Complementary Problem 3).

### Income shares and the Lorenz curve

Let's begin by working out the quantiles of the Pareto distribution.  The
starting point will be the CCDF:
\begin{eqnarray}
\overline{F}(x_q) & = & q\\
\left( \frac{x_q}{\xmin} \right)^{-\alpha+1} & = & q\\
\frac{x_q}{\xmin} & = & q^{1/(1-\alpha)}\\
x_q & = & \xmin q^{1/(1-\alpha)}
\end{eqnarray}

With this in hand, we can work out what fraction of total income goes to the
upper $q$ fraction of the population, call this $S(q)$:
\begin{eqnarray}
S(q) & = & \frac{\int_{x_q}^{\infty}{x \frac{\alpha-1}{\xmin}\left(\frac{x}{\xmin}\right)^{-\alpha} x}}{\Expect{X}}\\
& = & \frac{\frac{\alpha-1}\xmin^{\alpha-1}\int_{\xmin q^{1/(1-\alpha)}}^{\infty}{x^{-\alpha+1} dx}}{\xmin \frac{\alpha-1}{\alpha-2}}\\
& = & \xmin^{\alpha-2}\frac{\frac{1}{\alpha-2} \xmin^{-\alpha+2} q^{(\alpha-2)/(\alpha-1)}}{\frac{1}{\alpha-2}}\\
& = & q^{(\alpha-2)/(\alpha-1)}
\end{eqnarray}

Now by convention the Lorenz curve is plotted using _lower_ percentiles, say $p=1-q$, so we can say
\begin{equation}
s(p) = 1-(1-p)^{(\alpha-2)/(\alpha-1)}
\end{equation}

Here are some pictures:

```{r}
pareto.share <- function(p, alpha) { 1-(1-p)^((alpha-2)/(alpha-1)) }
par(pty="s")
curve(pareto.share(x, 3.5),
      xlim=c(0,1), ylim=c(0,1),
      xlab="Income percentile", ylab="Cumulative income share",
      lty=1)
curve(pareto.share(x, 2.5), add=TRUE, lty=2)
curve(pareto.share(x, 2.1), add=TRUE, lty=3)
legend("topleft", legend=c(3.5, 2.5, 2.1), lty=1:3, cex=0.5)
abline(a=0, b=1, col="grey")
```

We can see that as $\alpha$ decreases, the curve gets further and further away
from the diagonal, indicating increasing inequality.

(Note that this result will break down once $\alpha < 2$ --- why?)

### The Gini index

Recall, from Lecture 2, that the Gini index is determined by the area under
the Lorenz curve:
\begin{equation}
G = 1 - 2\times \text{area}
\end{equation}
The area under the Lorenz curve is straightforward to find from the income shares:
\begin{eqnarray}
\text{area under Lorenz curve} & = & \int_{0}^{1}{s(p) dp}\\
& = & \int_{0}^{1}{\left(  1-(1-p)^{(\alpha-2)/(\alpha-1)} \right) dp}\\
& = & 1 - \int_{0}^{1}{(1-p)^{(\alpha-2)/(\alpha-1)} dp}\\
& = & 1 - \int_{1}^{0}{ q^{(\alpha-2)/(\alpha-1)} (-1) dq}\\
& = & 1 - \int_{0}^{1}{ q^{(\alpha-2)/(\alpha-1)} dq}\\
& = & 1 - \frac{1}{\frac{\alpha-2}{\alpha-1} + 1}[1-0]\\
& = & 1 - \frac{1}{\frac{\alpha-2 + \alpha -1}{\alpha-1}}\\
& = & 1 - \frac{\alpha-1}{2\alpha -3}\\
& = & \frac{2\alpha -3 - \alpha +1}{2\alpha -3} = \frac{\alpha - 2}{2\alpha -3}
\end{eqnarray}
so
\begin{eqnarray}
G & = & 1 - 2  \frac{\alpha - 2}{2\alpha -3}\\
& =& \frac{2\alpha -3 - 2\alpha+4}{2\alpha -3} = \frac{1}{2\alpha-3}
\end{eqnarray}
This leads to the following relationship between $\alpha$ and the Gini
index:

```{r}
gini.from.alpha <- function(alpha) { 1/(2*alpha-3) }
curve(gini.from.alpha(x), from=2, to=10,
      xlab=expression(alpha),
      ylab="Gini index",
      main="Gini index of Pareto distributions")
```

## Lorenz curves and Gini indices for the log-normal distribution

Homework 1 quoted some results about the shape of the Lorenz curve and the
Gini index for the log-normal distribution, but without asking you to prove
them.  Here are quick sketches of the proofs.  (It's probably OK just
to skim these, at least for now.)

### Income shares and the Lorenz curve

Let's start by working out the (lower) quantiles of the $\mathcal{LN}(\mu, \sigma^2)$ distribution:
\begin{eqnarray}
\Prob{X \leq x_p} & = & p\\
\Prob{X \leq \log{x_p}} & = & p\\
\Phi((\log{x_p} - \mu)/\sigma)) & = & p\\
\log{x_p} & = & \mu + \sigma \Phi^{-1}(p)\\
x_p & = & e^{\mu} \myexp{\sigma \Phi^{-1}(p)}
\end{eqnarray}

Now what's the share of income going to the lower $p$ fraction of the population?
\begin{eqnarray}
s(p) & = & \frac{\int_{0}^{x_p}{x f(x) dx}}{\Expect{X}}\\
& = & \frac{\int_{0}^{e^{\mu} \myexp{\sigma \Phi^{-1}(p)}}{\frac{1}{\sqrt{2\pi\sigma^2}} \myexp{-\frac{1}{2}\frac{(\mu - \log{x})^2}{\sigma^2} dx}}}{e^\mu e^{\sigma^2/2}}\\
\end{eqnarray}
Set $z = \log{x}$, so $dz/dx = 1/x$ but $dx = x dz = e^{z} dz$; make $z$ the new variable of integration and change the limits accordingly:
\begin{eqnarray}
\lefteqn{\int_{-\infty}^{e^{\mu} \myexp{\sigma \Phi^{-1}(p)}}{\frac{1}{\sqrt{2\pi\sigma^2}} \myexp{-\frac{1}{2}\frac{(\mu - \log{x})^2}{\sigma^2}} dx} } & & \\
& = & \int_{-\infty}^{\mu + \sigma \Phi^{-1}(p)}{ \frac{1}{\sqrt{2 \pi \sigma^2}} \myexp{-\frac{1}{2}\frac{(\mu - z)^2}{\sigma^2}} e^{z} dz}\\
& = & \int_{-\infty}^{\mu + \sigma\Phi^{-1}(p)}{\frac{1}{\sqrt{2\pi\sigma^2}} \myexp{-\frac{1}{2\sigma^2}(\mu^2 - 2\mu z - 2 z \sigma^2 + z^2)} dz}\\
& = & \int_{-\infty}^{\mu + \sigma\Phi^{-1}(p)}{\frac{1}{\sqrt{2\pi\sigma^2}} \myexp{-\frac{1}{2\sigma^2}(\mu^2 - 2(\mu +\sigma^2) z + z^2)} dz}\\
& = & \int_{-\infty}^{\mu + \sigma\Phi^{-1}(p)}{\frac{1}{\sqrt{2\pi\sigma^2}} \myexp{-\frac{1}{2\sigma^2}((\mu+\sigma^2)^2 - 2(\mu + \sigma^2) z + z^2 -\sigma^4 - 2\mu\sigma^2)} dz}\\
& = & \int_{-\infty}^{\mu + \sigma\Phi^{-1}(p)}{\frac{1}{\sqrt{2\pi\sigma^2}} \myexp{\sigma^2/2 + \mu}\myexp{-\frac{1}{2\sigma^2}((\mu+\sigma^2 -z)^2} dz}\\
& = & e^{\sigma^2/2 + \mu}\int_{-\infty}^{\mu + \sigma\Phi^{-1}(p)}{\frac{1}{\sqrt{2\pi\sigma^2}} \myexp{-\frac{1}{2\sigma^2}((\mu+\sigma^2 -z)^2} dz}
\end{eqnarray}
At the end of all this manipulation ("completing the square"), what we're left
with as the integrand is the pdf of a Gaussian with mean $\mu+\sigma^2$ and
variance $\sigma^2$, so this becomes
\begin{eqnarray}
\lefteqn{e^{\sigma^2/2 + \mu}\Phi\left( \frac{(\mu + \sigma\Phi^{-1}(p))-(\mu+\sigma^2)}{\sigma}\right)} & & \\
& = & e^{\sigma^2/2+\mu}\Phi(\Phi^{-1}(p) - \sigma)
\end{eqnarray}
Thus
\begin{eqnarray}
s(p) & = & \frac{e^{\sigma^2/2+\mu}\Phi(\Phi^{-1}(p) - \sigma)}{e^{\mu + \sigma^2/2}}\\
& = & \Phi(\Phi^{-1}(p) - \sigma)
\end{eqnarray}
as claimed (without proof) in Homework 1.

```{r}
lorenz.lnorm <- function(p, sigma) {
    pnorm(qnorm(p) - sigma)
}

par(pty="s")
curve(lorenz.lnorm(x, sigma=1),
      from=0, to=1, ylim=c(0,1),
      xlab="Income percentile", ylab="Cumulative income share", lty=1)
curve(lorenz.lnorm(x, sigma=sqrt(10)), add=TRUE, lty=2)
curve(lorenz.lnorm(x, sigma=sqrt(.1)), add=TRUE, lty=3)
legend("topleft", legend=c(expression(sigma^2==1),
                           expression(sigma^2==10),
                           expression(sigma^2==0.1)),
       lty=1:3, cex=0.5)
abline(0, 1, col="grey")
```

Notice that $\mu$ doesn't appear in $s(p)$ at all, so two log-normal
distributions with different $\mu$s but the same $\sigma$ will have the same
Lorenz curve.  This makes some sense: changing $\mu$ will change everyone's
income by the same factor, and that's precisely the kind of change which the
Lorenz curve ignores.

### Gini index

This $2\Phi(\sigma/\sqrt{2})-1$, but you're tired of reading me do integrals by
now, so I'll make this Complementary Problem 3.  Notice that since the Lorenz
curve for a log-normal only involves $\sigma$ and not $\mu$, the Gini index
couldn't involve $\mu$.


```{r}
curve(2*pnorm(sqrt(x)/sqrt(2))-1, from=0, to=10,
      xlab=expression(sigma^2),
      ylab="Gini index",
      main="Gini index of log-normal distributions")
```

\clearpage


# Complementary Problems

These are to think through or practice on, not to hand in.

1. The code for this lecture downloaded income data for multiple countries and
   years (as in Lecture 2), but I've only shown the pdf of the income
   distribution for the US in 2019 (as in Lecture 1).  Pick another country and
   year, say Finland in 2005, and make similar plots of the pdf of income.  Do
   this again for a bunch of other countries and years.
2. _Some calculus with the log-normal distribution_
    a. Show that the mode of the log-normal distribution is $e^{\mu-\sigma^2}$.
	    There are a couple of ways to do this, but the must straightforward is
		to solve $\frac{df}{dx} = 0$ for $x$.
	b. Show that $\Expect{X} = e^{\mu + \sigma^/2}$.  Again, there are a couple
	   of ways to do this, including brute-force integration.
3. Show that the Gini index of a log-normal distribution is $2\Phi(\sigma/\sqrt{2})-1$.  One way to do this is to find the area under the Lorenz curve.



