Talent vs Luck simulation

tvl.R

library(checkpoint)
checkpoint("2018-02-25")
library(ggplot2)

# number of people
N <- 1000
# probability of event interception
P_E <- 0.075
# probability of lucky event
P_L <- 0.5
# initial capital
C_0 <- 10.
# total time steps
T_ <- 80
# talent mean
M_T <- 0.6
# talent standard deviation
SD_T <- 0.1

talents <- rnorm(N, M_T, SD_T)
c_all <- matrix(0, T_+1, N)
c_all[1, ] <- rep(C_0, times=N)
event_all <- matrix(0, T_, N)
for(i in seq(T_)){
    # whether an event happens
    event_happens <- rbinom(N, 1, P_E)
    # whether an event is lucky/good
    event_is_good <- rbinom(N, 1, P_L)
    # whether a person can profit from a lucky event
    event_good_profited <- rbinom(N, 1, talents)
    # The effect of the hypothetical event
    event_effect_hypo <- ((event_is_good & event_good_profited) * 3 + 1) / 2
    # The actual effect at this time step
    event_effect_actual <- event_effect_hypo * event_happens
    event_effect_actual[event_effect_actual==0] <- 1
    event_all[i, ] <- event_effect_actual
    c_all[(i+1),] <- c_all[i,] * event_effect_actual
} 

ggplot(data.frame(talents=talents), aes(x=talents)) + geom_histogram(binwidth=0.02) + 
  geom_vline(xintercept=0.6, linetype=2) + geom_vline(xintercept=0.7, linetype=3) + 
  ylab("number of individuals") +
  geom_vline(xintercept=0.5, linetype=3) + theme_bw() + theme(text=element_text(size=14))

# bins <- cut(c_all[T_+1,], breaks=1000, include.lowest = T, right=FALSE)

c_hist <- hist(c_all[T_+1,], breaks=100, plot=F)
plot(
  (c_hist$breaks[1:(length(c_hist$breaks)-1)] + c_hist$breaks[2:length(c_hist$breaks)])/2,
  c_hist$counts, log="y", type='h', lwd=2, lend=2, 
  xlab="capital/success", 
  ylab="number of individuals")


min(log10(c_all[T_+1,]))
max(log10(c_all[T_+1,]))
c_hist <- hist(pmax(log10(c_all[T_+1,]),-2), 
               breaks=c(
                 floor(min(log10(c_all[T_+1,]))), 
                 seq(1, ceiling(max(log10(c_all[T_+1,]))), 0.5)), plot=F)
ggplot(data.frame(
    x=seq(1, ceiling(max(log10(c_all[T_+1,]))), 0.5), cnt=log10(c_hist$count)
  ), aes(x=x, y=cnt)) + 
  geom_point() +
  geom_smooth(method = "lm", se = T, col="red") + 
  scale_x_continuous("capital/success", seq(1, ceiling(max(log10(c_all[T_+1,]))), 1), 
  labels=10 ^ seq(1, ceiling(max(log10(c_all[T_+1,]))), 1)) +
  scale_y_continuous("number of individuals", seq(0, 3, 1), 
                     labels=10 ^ seq(0, 3, 1)) +  
  theme_bw() + theme(text=element_text(size=14))

ggplot(data.frame(x=c_all[T_+1,], y=talents), aes(x=x, y=y)) + 
  geom_point(alpha=0.25) + ylab("talent") + xlab("capital/success") + ylim(0, 1) +
  scale_x_log10(breaks=c(1, 10, 100, 1000)) + theme_bw() + 
  theme(text=element_text(size=14))
  
ggplot(data.frame(y=c_all[T_+1,], x=talents), aes(x=x, y=y)) + 
  geom_point(alpha=0.5) + xlab("talent") + ylab("capital/success") + theme_bw() + 
  theme(text=element_text(size=14))

ggplot(data.frame(y=c_all[T_+1,], x=talents), aes(x=x, y=y)) + 
  geom_point(alpha=0.5) + xlab("talent") + ylab("capital/success") +
  scale_y_log10(breaks=c(1, 10, 100)) + theme_bw() + 
  theme(text=element_text(size=14))

c_final_sorted <- sort(c_all[T_,])
sum(c_final_sorted)
sum(c_final_sorted[1:(N*0.8)]) / sum(c_final_sorted)
sum(c_final_sorted[(N*0.8):N]) / sum(c_final_sorted)

# Performance of the Talented 
talented <- c_all[T_,][talents > 0.7]
length(talented)
sum(talented > 10) / length(talented)
sum(c_all[T_,] > 10) / N