# Given summary statistics
mean_income <- 84000
median_income <- 62000
# Calculate lognormal parameters
mu <- log(median_income)
sigma <- sqrt(2 * log(mean_income / median_income))
# Simulate incomes
set.seed(123)
n <- 10000
income <- rlnorm(n, meanlog = mu, sdlog = sigma)
# Plot histogram (density scale)
hist(income,
probability = TRUE,
breaks = 50,
col = "lightgray",
border = "white",
main = "Simulated Lognormal Income Distribution",
xlab = "Income ($)",
ylab = "Density",
xaxt = "n")
# Add x-axis with comma formatting
axis(1,
at = pretty(income),
labels = format(pretty(income),
big.mark = ",",
scientific = FALSE))
# Overlay the theoretical lognormal density in red
curve(dlnorm(x, meanlog = mu, sdlog = sigma),
from = 0,
to = max(income),
add = TRUE,
col = "red",
lwd = 3)
# Verify simulation
cat("Sample Mean: $", format(round(mean(income), 0),
big.mark = ","), "\n")
## Sample Mean: $ 83,856
cat("Sample Median: $", format(round(median(income), 0),
big.mark = ","), "\n")
## Sample Median: $ 61,466
print("thr")
## [1] "thr"
# Probability of earning more than $100,000
1 - plnorm(100000, meanlog = mu, sdlog = sigma)
## [1] 0.2698098
print("sim")
## [1] "sim"
mean(income > 100000)
## [1] 0.2718
print("thr")
## [1] "thr"
# Probability of earning more than $1,000,000
1 - plnorm(1000000, meanlog = mu, sdlog = sigma)
## [1] 0.0001799027
print("sim")
## [1] "sim"
mean(income > 1000000)
## [1] 2e-04
amt_yr <- 32000
mean(income < 32000)
## [1] 0.1977
plnorm(32000, meanlog = mu, sdlog = sigma)
## [1] 0.1980327
# Federal tax brackets (single filer)
brackets <- c(0, 12400, 50400, 105700, 201775, 256225, 640600, Inf)
labels <- c("10%", "12%", "22%", "24%",
"32%", "35%", "37%")
tax_bracket <- cut(income,
breaks = brackets,
labels = labels,
right = TRUE,
include.lowest = TRUE)
# Number of people in each bracket
table(tax_bracket)
## tax_bracket
## 10% 12% 22% 24% 32% 35% 37%
## 189 3747 3593 1835 294 327 15
# Percentage of people in each bracket
round(prop.table(table(tax_bracket)) * 100, 2)
## tax_bracket
## 10% 12% 22% 24% 32% 35% 37%
## 1.89 37.47 35.93 18.35 2.94 3.27 0.15
# Marginal tax rates
rates <- c(0.10, 0.12, 0.22, 0.24, 0.32, 0.35, 0.37)
# Function to compute federal tax
calc_tax <- function(x, brackets, rates) {
tax <- 0
for(i in 1:length(rates)) {
lower <- brackets[i]
upper <- brackets[i + 1]
if(x > lower) {
taxable <- min(x, upper) - lower
tax <- tax + taxable * rates[i]
} else {
break
}
}
return(tax)
}
# Apply to every simulated income
amt_yr <- sapply(income, calc_tax,
brackets = brackets,
rates = rates)
# Average annual tax
mean(amt_yr)
## [1] 14919.69
# Average effective tax rate
mean(amt_yr / income)
## [1] 0.144844
# First few observations
head(data.frame(
Income = round(income, 0),
Tax = round(amt_yr, 0),
EffectiveRate = round(100 * amt_yr / income, 2)
))
## Income Tax EffectiveRate
## 1 40058 4559 11.38
## 2 51819 6112 11.80
## 3 208908 43306 20.73
## 4 65502 9122 13.93
## 5 68573 9798 14.29
## 6 235980 51970 22.02
after_tax_income <- income - amt_yr
# Histogram of after-tax income
hist(after_tax_income,
probability = TRUE,
breaks = 50,
col = "lightblue",
border = "white",
main = "After-Tax Income Distribution",
xlab = "After-Tax Income ($)",
ylab = "Density",
xaxt = "n")
# Format x-axis with commas
axis(1,
at = pretty(after_tax_income),
labels = format(pretty(after_tax_income),
big.mark = ",",
scientific = FALSE))
# Overlay kernel density estimate
lines(density(after_tax_income),
col = "red",
lwd = 3)
after_tax_mn <- mean(after_tax_income)
after_tax_md <- median(after_tax_income)
mu <- log(after_tax_mn)
sigma <- sqrt(2 * log(after_tax_mn / after_tax_md))
print("thr")
## [1] "thr"
# Probability of earning more than $100,000
1 - plnorm(100000, meanlog = mu, sdlog = sigma)
## [1] 0.302465
print("sim")
## [1] "sim"
mean(income > 100000)
## [1] 0.2718
print("thr")
## [1] "thr"
# Probability of earning more than $1,000,000
pct <- 1 - plnorm(1000000, meanlog = mu, sdlog = sigma)
print("sim")
## [1] "sim"
mean(income > 1000000)
## [1] 2e-04
pct *100
## [1] 0.009979663