food things and veg

food.r

library(R.oo)
library(plyr)

# Main function. Given the number of new vegetarians (in millions),
# returns the number of people brought into or out of poverty.
changeInPoor = function(milVeg) {
  # Demand in the F&P paper was 61 million bushels
  # We find the % change in feed demand, and use that with the elasticity
  # from F&P to find % change in price
  relDemandChange = vegChange(milVeg) / 61
  relPriceChange = .0868 * relDemandChange # dPrice / dSupply = .0868
  #relPriceChange = absPriceChange / 7.36  # Corn price per bushel as of 2/3/13
  relPriceChange = relPriceChange / 3     # Corn + Soy = 1/3 of staples
  print(paste("Change in feed demand: ", relDemandChange * 100, "%", sep=''))
  print(paste("Change in price: ", relPriceChange*100,"%", sep=""))
  poorChange  = newPoor(1 + relPriceChange)
  print(paste("Change in number of poor:", poorChange))
  
  # Goklany estimate on DALY averted by bringing people out of poverty
  # from http://www.jpands.org/vol16no1/goklany.pdf
  dalysAverted = (-183000 / 1000000) * poorChange
  print(paste("DALYs averted:", dalysAverted))
  
  print(paste("Cost per DALY:", (milVeg * 1000000 * 11) / dalysAverted, sep=" $"))
}

# Given a change in food price, calculates the change
# in number of poor (< $2/day)
# E.g. newPoor(.99) finds the number of poor if food prices dropped 1%
newPoor = function(foodChange, gidd = loadGidd()) {
  oldPoor = numPoor(gidd)
  cvs = mapply(food, foodChange, gidd$foodshare, .2) 
  # e.g. if compensating variation is 1.2 times your income,
  # we say real income falls to .8 of nominal
  incomeChanges = (2 - cvs)
  gidd$consincPPP05 = gidd$consincPPP05 * incomeChanges
  newPoor = numPoor(gidd)
  return(newPoor - oldPoor)
}

# Fraction of population making less than $37.5/mo.
# Encapsulated in case we have a more complex eqn in the future
numPoor = function(gidd = loadGidd()) {
  # return(sum(gidd$pop[gidd$consincPPP05 <= 37.5]))
  poorPerCountry = ddply(gidd, .(country),
        function(cntry){
          totInPov = 0
          for(i in 1:nrow(cntry)){
            # Assume a uniform distribution of income between endpoints in 
            # a vintile
            #
            # Assume that the lowest income is 1/2 of mean and the highest income is
            # twice mean
            if (i == 1) {lower = 0} 
            else {lower = .5 * cntry[i-1,"consincPPP05"]}
            if (i == nrow(cntry)) {upper = 2 * cntry[i,"consincPPP05"]}
            else {upper = cntry[i+1, "consincPPP05"]}
            
            # Find the fraction of people in this vintile below the poverty
            # line of 37.5
            if(lower > 37.5) next # No one in the vintile makes <= 37.5
            fracUnderPov = (37.5 - lower) / (upper - lower)
            totInPov = totInPov + fracUnderPov * cntry[i, "pop"]
          }
          return(totInPov)
        })
  return(sum(poorPerCountry$V1))
}

# Dessus et al. equation for change in real income as a function
# of food price changes.
food = function(priceChange, foodCons, priceElast) {
  lpf = log(priceChange)
  lpn = log(1)  # Assume non-food stays constant
  
  wFood = foodCons
  wNonFood = 1 - wFood
  eFF = priceElast
  eFN = priceElast
  eNN = priceElast
  lnC = wFood * lpf + wNonFood * lpn +
           .5 * (wFood * eFF * lpf * lpf +
                   wFood * eFN * lpf * lpn +
                   wFood * eFN * lpf * lpn +
                   wFood * eNN * lpn * lpn)
  
   #lnC = wFood * lpf + .5 * (wFood * eFF * lpf * lpf)
  return(exp(lnC))
}

# From http://www.aae.wisc.edu/renk/library/Effect%20of%20Ethanol%20on%20Corn%20Price.pdf
feedDemandChange = function(broilers, hogs, cows) {
  Const = -1.4863
  PCt = -.2978
  Psmt = .1813
  Bt = .1999
  COFt = .3885
  Ht = .4423
  D1 = .3851
  D2 = .2410
  D3 = .1212
  
  return(Bt * broilers + Ht * hogs + COFt * cows)
}

cornPriceChange = function(broilers, hogs, cows) {
  relDemandChange = (feedDemandChange(broilers, hogs, cows)) / 61
  return(.0868 * relDemandChange)
}

# Change in corn demand for feed, based on the Counting Animals estimate of
# animals saved
vegChange = function(milVeg) {
  return(feedDemandChange(-28.637 * milVeg, -.44 * milVeg, -.132 * milVeg))
}

# Helper function to parse tableBig
formatTable = function(tableStr){
  tableStr = trim(tableStr)
  tableStr = gsub('\n[ ]+','\n', tableStr)
  tableStr = gsub("[ ]+",',', tableStr)
  table = read.csv(textConnection(tableStr), header = F)
}

loadGidd = function() {
  gidd = read.csv('C:\\Users\\ben\\Downloads\\GlobalDist - global_dist.csv')
  # For rows that are missing the share of income going to food, set it to be 
  # the average of other nations with similar incomes
  foodShares = ddply(gidd[is.na(gidd$foodshare),], 
        # country & vintile uniquely identify a row
        .(country, vintile), 
        function(row){
          similarRows = gidd[gidd$consincPPP05 > 0.8 * row$consincPPP05 &
                             gidd$consincPPP05 < 1.2 * row$consincPPP05 &
                             !is.na(gidd$foodshare),
                       ]
          return(c(mean(similarRows$foodshare), nrow(similarRows)))
  })
  gidd$foodshare[is.na(gidd$foodshare)] = foodShares$V1
  return(gidd)
}