Create a gist now

Instantly share code, notes, and snippets.

@Xodarap /food.r
Last active Dec 12, 2015

What would you like to do?
food things and veg
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)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment