Simulations

Published:

Entrepreneurship Model Simulation

Cost Function:

\[ c(t) = \begin{cases} dt \text{ if } t \leq T \\ Dt \text{ if } t > t \end{cases} \]

  • d = per period cost on leave \(\forall t \leq T\)
  • (D-d) = discrete increase \(t > T\)
  • J = jump in cost

Benefit Function

\[ b(t) = st ^{1/\alpha}\gamma \]

Where \(\gamma = \Gamma(1 - \frac{1}{\alpha})\)

R gamma function returns \[ \Gamma(x)= \int_{0}^{\infty} t^{x-1} e^{-t} dt \]

Utility function:

\[ u(t) = st ^{1/\alpha}\gamma - c(t) \]

R function to plot Cost and Benefit

  • and all relevant values in the subtitle + caption
plot_functions <- function(d,D,FinalT,s,alpha,
                           plotUtility=FALSE) {
  # casewise cost function
  cost <- function(t){
    ifelse(t <= FinalT, d*t,
    ifelse(t > FinalT,  D*t, NA ))
  }
  gamma = gamma(1-(1/alpha))
  # benefit function
  benefit <- function(t) {s*t^(1/alpha)*gamma}
  # guess that cost and benefit intersect between t = 0.001 and
  # 100; ignore trivial solution at t = 0
  intersection = uniroot(
    function(t) benefit(t)-cost(t),c(0.001,100),
    extendInt='yes')$root
  # critical values for solution regimes
  st_alpha_gamma = s*FinalT^(1/alpha)*gamma
  dT = d*FinalT
  DT = D*FinalT
  # utility function

  utility <- function(t) {benefit(t) - cost(t)}
  t_star = optimise(utility, interval=c(0.001, 100),
                    maximum=TRUE)$maximum

  # optima
  crit_1    = ((s*gamma)/d)^(alpha/(alpha-1))
  optimum_1 = ((s*gamma)/(d*alpha))^(alpha/(alpha-1))

  crit_2    = ((s*gamma)/D)^(alpha/(alpha-1))
  optimum_2 = ((s*gamma)/(D*alpha))^(alpha/(alpha-1))

  # ifelse (condition, value if true,
    # condition if false(condition 2, value if true, value if false))
  t_optimum = ifelse(crit_1<FinalT, optimum_1,
                ifelse(crit_2>FinalT,optimum_2, FinalT))

  if (plotUtility==FALSE){
    fns <- ggplot(data.frame(t=c(0,18)), aes(t)) +
    stat_function(fun=cost, aes(colour= 'cost'),
                  geom='step', n = 1000)+
    stat_function(fun=benefit, aes(colour = 'benefit'))
  }
  else {
    fns <- ggplot(data.frame(t=c(0,18)), aes(t)) +
    stat_function(fun=cost, aes(colour= 'cost'),
                  geom='step', n = 1000)+
    stat_function(fun=benefit, aes(colour = 'benefit'))+
    stat_function(fun=utility, aes(colour = 'utility'))
  }
  plot = fns + labs(title = "Cost and Benefit Functions",
       subtitle = paste0(
        'T=', FinalT,
        ';d=', d,
        ';D=',D,
        ';s=',s,
        # ';gamma=',format(gamma,digits=2),
        ';alpha=',format(alpha,digits=2)),
       x='Time',
       caption=paste0('intersection=',
                      format(intersection,digits=3),
                      ';t*=',format(t_star,digits=3)
                      # '\n ;(sg/d)^r=',format(crit_1,digits=3),
                      # ';(sg/da)^r=',format(optimum_1,digits=3),
                      # '\n ;(sg/D)^r=',format(crit_2,digits=3),
                      # ';(sg/Da)^r=',format(optimum_2,digits=3)
                      # ';numerical t*=',
                      # format(t_star,digits=3)
                      ))+
      scale_colour_brewer(palette='Set1')+
      labs(colour="Function")+
      scale_x_continuous(breaks=c(0,3,6,9,12,15,18))+
      geom_vline(aes(xintercept=t_star),
                color="black", linetype="dashed")
}

Plots:

Comparative Statics : Parameter changes within same solution regime (\(sT^{1/\alpha}\gamma < dT\))

baseline

plot = plot_functions(d=1, D=1.2,
              FinalT = 12, s = 0.5, alpha = 1.5 )
print(plot)

Increasing cost of experimenting

Increasing d

plot= plot_functions(d = 1.1, D = 1.2,
              FinalT = 12, s = 0.5, alpha = 1.5 )
print(plot)

Increasing value of experimenting

Increasing \(\gamma\) (decreasing \(\alpha\) )

plot= plot_functions(d=1, D=1.2,
              FinalT = 12, s = 1, alpha = 1.47 )
print(plot)

Increasing s

plot= plot_functions(d=1, D=1.2,
              FinalT = 12, s = 1.5, alpha = 1.5 )
print(plot)

Comparing across solution regimes

I now increase s s.t. \(dT < sT^{1/\alpha} < DT+J\)

plot= plot_functions(d=1, D=1.2,
              FinalT = 12, s = 2.4,alpha = 1.5)
print(plot)

Holding cost constant and adding multiple benefit functions

d = 1
D = 2
FinalT = 12

alpha = 1.7
gamma = gamma(1-(1/alpha))
# benefit functions
benefit1 <- function(t) {1.35*t^(1/alpha)*gamma}
benefit2 <- function(t) {1.5*t^(1/alpha)*gamma}
benefit3 <- function(t) {1.65*t^(1/alpha)*gamma}
benefit4 <- function(t) {2*t^(1/alpha)*gamma}


plot= plot_functions(d = 1, D = 1.2,
              FinalT = 12, s = 1,alpha = 1.7)+
      stat_function(fun=benefit1, aes(colour = 's=1.35'))+
      stat_function(fun=benefit2, aes(colour = 's=1.5'))+
      stat_function(fun=benefit3, aes(colour = 's=1.65'))+
      stat_function(fun=benefit4, aes(colour = 's=2'))+
        ylim(0,35)
print(plot)

Spike at t = T

Draw s from a uniform distribution, store optimal t, repeat

leave_time_dist <-function(d,D,FinalT,alpha,n,ub,lb){
  # draw n values of s from a uniform distribution between 1 and 15
  u = runif(n,lb,ub)
  # initialise empty array of solutions
  t_stars = numeric()
  t_stars_a = numeric()
  for (s in u){
    # cost function
    cost <- function(t){
      ifelse(t <= FinalT, d*t,
      ifelse(t > FinalT,  D*t, NA ))
    }
    gamma = gamma(1-(1/alpha))
    # benefit function
    benefit <- function(t) s*t^(1/alpha)*gamma
    # utility function = difference between the two
    utility <- function(t) benefit(t) - cost(t)
    # find utility maximiser value of t
    soln = optimise(utility, interval=c(0.001, 150),
                    maximum=TRUE)$maximum
    # append to array
    t_stars = c(t_stars,soln)

    # optima
    crit_1    = ((s*gamma)/d)^(alpha/(alpha-1))
    optimum_1 = ((s*gamma)/(d*alpha))^(alpha/(alpha-1))

    crit_2    = ((s*gamma)/D)^(alpha/(alpha-1))
    optimum_2 = ((s*gamma)/(D*alpha))^(alpha/(alpha-1))

    # ifelse (condition, value if true,
    # condition if false(condition 2, value if true,
    # value if false))

    t_optimum = ifelse(crit_1<FinalT, optimum_1,
                  ifelse(crit_2>FinalT,optimum_2, FinalT))

    t_stars_a = c(t_stars_a,t_optimum)

  }
  # store as dataframe, pass to ggplot
  dist = cbind.data.frame(u,t_stars,t_stars_a)
  names(dist)
  hist_u = ggplot(data = dist,aes(u)) +
    geom_freqpoly(binwidth=0.05) +
    labs(title='Distribution of s',
         subtitle=paste0(n,' draws from [',
                         lb,',',ub,']'),
         x = 'value of experimenting (s)')

  hist_t = ggplot(data = dist,aes(t_stars)) +
    geom_freqpoly(binwidth=0.05)+
    labs(title='Distribution of t*',
         x = 'leave time taken (t*)')+
    scale_x_continuous(breaks=c(3,6,9,12,15,18))+
    theme_bw() +
    coord_cartesian(ylim=c(0,n*0.01),
                    xlim=c(9,15))

  multiplot(hist_u,hist_t)

  }
leave_time_dist(d = 1, D = 1.2,
              FinalT = 12,alpha = 1.6,
              n = 100000, lb = 1, ub = 7)

Probability to choose Entrepreneurship

\[ Pr (\pi > w_{t*}) = 1 - \exp [-t^* (\frac{1}{s} (1-\delta)^{t*\mathbb{1}_{t^*>T}} w_1)^{-\alpha}] \]

prob_choose_entrepreneurship <- function(FinalT,alpha,n,
                                         w1,r=0.95,delta=0.25,theta=0.2,
                                         ub,lb){
  d = r + theta * delta # 0.95 + 0.05 = 1 by default
  D = r + delta         # 0.95 + 0.25 = 1.2 by default

  u = runif(n,lb,ub)
  # initialise empty array of solutions
  t_stars = numeric()
  prob = numeric()
  for (s in u){
    # cost function
    cost <- function(t){
      ifelse(t <= FinalT, d*t,
      ifelse(t >  FinalT, D*t, NA ))
    }
    gamma = gamma(1-(1/alpha))
    # benefit function
    benefit <- function(t) s*t^(1/alpha)*gamma
    # utility function = difference between the two
    utility <- function(t) benefit(t) - cost(t)
    # find utility maximiser value of t
    t_star = optimise(utility, interval=c(0.001, 150),
                    maximum=TRUE)$maximum
    # append to array
    t_stars = c(t_stars,t_star)

    indicator <- ifelse(t_star <= FinalT, 0,
                        ifelse(t_star > FinalT,  1, NA ))

    pr = 1 - exp(-t_star*((1/s)*(1-delta)^(t_star*indicator)*w1)^(-alpha))
    prob = c(prob,pr)

  }
  # store as dataframe, pass to ggplot
  dist = cbind.data.frame(u,t_stars,prob)

  hist_u = ggplot(data = dist,aes(u)) +
    geom_freqpoly(binwidth=0.05) +
    labs(title='Distribution of s',
         subtitle=paste0(n,' draws from [',
                         lb,',',ub,']'),
         x = 'value of experimenting (s)')

  hist_prob = ggplot(data = dist,aes(prob)) +
    geom_histogram(binwidth=0.01)+
    labs(title='Probability of choosing Entrepreneurship',
         x = 'Probability')+
    theme_bw()

  dist %>% filter(t_stars<=FinalT) %>%
    rename(t_star = t_stars) ->
    leq_finalT

  scatterplot = ggplot(data =leq_finalT,aes(t_star,prob))+
    geom_point()+
    scale_x_continuous(breaks=c(0,3,6,9,12))+
    labs(x = 't*', y = 'prob')

  multiplot(hist_u,hist_prob,scatterplot)
}
prob_choose_entrepreneurship(FinalT = 12,alpha = 1.6,n = 100000,
              w1 = 5,
              lb = 0.01, ub = 4)