<jupyter_start><jupyter_text>Minimum Cost Diet<jupyter_code>import ndb
import pandas as pd
import warnings

DFs = []
#for csv in CSVs:  # Uncomment to use a list of csv files as inputs
#  DFs.append(pd.read_csv(csv,dtype=str))

try:
    if len(SHEETs):
        for ID, RANGE_NAME in SHEETs:
            try:
                if "docs.google.com" in ID:
                    sheet = "%s&output=csv" % ID
                else:
                    sheet = "https://docs.google.com/spreadsheet/ccc?key=%s&output=csv" % ID
                DFs.append(pd.read_csv(sheet))
            except ParserError:
                warnings.warn("Can't read sheet at https://docs.google.com/spreadsheets/d/%s.\nCheck Sharing settings, so that anyone with link can view?" % ID)
except NameError: # SHEETS not defined?
    pass

df = pd.concat(DFs,ignore_index=True,sort=False)
# Some columns which ought to be numeric are actually str; convert them
df['Price'] = df['Price'].astype(float)
df['Quantity'] = df['Quantity'].astype(float)
df["Units"] = df["Units"].astype(str)
df["Value"] = df["Value"].dropna().astype(int).astype(str)
df["NDB"][df['Food'].str.contains("Milk shakes and other dairy drinks")] = str("01110")
df["NDB"] = df["NDB"].dropna()

df
values = df["Value"].tolist()
names = df["Food"].tolist()
value_name_dict = dict(zip(values,names))
user = "atwater"
apikey = {'atwater':"<KEY>"}
age_list = ["C 1-3", "F 4-8", "M 4-8", "F 9-13", "M 9-13", "F 14-18", "M 14-18", "F 19-30", "M 19-30",\
            "F 31-50", "M 31-50", "F 51+", "M 51+"]
#http://www.censusscope.org/us/chart_age.html
age_breakdown = [6.82, 3.56, 3.74, 3.56, 3.74, 3.49, 3.69, 6.71, 6.92,\
                15.32, 15.15, 15.01, 12.32]
    
    
import ndb
D = {}

for food in  df.Food.tolist():
    try:
        NDB = df.loc[df.Food==food,:].NDB
        D[food] = ndb.ndb_report(apikey[user],NDB).Quantity
    except AttributeError: 
        warnings.warn("Couldn't find NDB Code %s for food %s." % (food,NDB))
        

D = pd.DataFrame(D,dtype=float)

D = D.reindex(sorted(D.columns), axis=1)

D
# Convert food quantities to NDB units
df['NDB Quantity'] = df[['Quantity','Units']].T.apply(lambda x : ndb.ndb_units(x['Quantity'],x['Units']))

# Now may want to filter df by time or place--need to get a unique set of food names.
df['NDB Price'] = df['Price']/df['NDB Quantity']

df.dropna(how='any') # Drop food with any missing data

# To use minimum price observed
Prices = df.groupby('Food')['NDB Price'].min()
# Choose sex/age group:
group = "M 19-30"

# Define *minimums*
bmin = pd.read_csv('./diet_minimums.csv').set_index('Nutrition')[group]

# Define *maximums*
bmax = pd.read_csv('./diet_maximums.csv').set_index('Nutrition')[group]

from  scipy.optimize import linprog as lp
import numpy as np

tol = 1e-6 # Numbers in solution smaller than this (in absolute value) treated as zeros

c = Prices.apply(lambda x:x.magnitude).dropna()

# Compile list that we have both prices and nutritional info for; drop if either missing
use = list(set(c.index.tolist()).intersection(D.columns.tolist()))
c = c[use]

# Drop nutritional information for foods we don't know the price of,
# and replace missing nutrients with zeros.
Aall = D[c.index].fillna(0)

# Drop rows of A that we don't have constraints for.
Amin = Aall.loc[bmin.index]

Amax = Aall.loc[bmax.index]

# Minimum requirements involve multiplying constraint by -1 to make <=.
A = pd.concat([-Amin,Amax])

b = pd.concat([-bmin,bmax]) # Note sign change for min constraints

# Now solve problem!
result = lp(c, A, b, method='interior-point')

# Put back into nice series
diet = pd.Series(result.x,index=c.index)

print("Cost of diet for %s is $%4.2f per day." % (group,result.fun))
print("\nYou'll be eating (in 100s of grams or milliliters):")
print(diet[diet >= tol])  # Drop items with quantities less than precision of calculation.

tab = pd.DataFrame({"Outcome":np.abs(A).dot(diet),"Recommendation":np.abs(b)})
print("\nWith the following nutritional outcomes of interest:")
print(tab)

print("\nConstraining nutrients are:")
excess = tab.diff(axis=1).iloc[:,1]
print(excess.loc[np.abs(excess) < tol].index.tolist())
age_list = ["C 1-3", "F 4-8", "M 4-8", "F 9-13", "M 9-13", "F 14-18", "M 14-18", "F 19-30", "M 19-30",\
            "F 31-50", "M 31-50", "F 51+", "M 51+"]
lowest_price_diet = []
for age in age_list:
    group = age
    bmin = pd.read_csv('./diet_minimums.csv').set_index('Nutrition')[group]
    bmax = pd.read_csv('./diet_maximums.csv').set_index('Nutrition')[group]
    tol = 1e-6
    c = Prices.apply(lambda x:x.magnitude).dropna()
    use = list(set(c.index.tolist()).intersection(D.columns.tolist()))
    c = c[use]
    Aall = D[c.index].fillna(0)
    Amin = Aall.loc[bmin.index]
    Amax = Aall.loc[bmax.index]
    A = pd.concat([-Amin,Amax])
    b = pd.concat([-bmin,bmax])
    result = lp(c, A, b, method='interior-point')
    diet = pd.Series(result.x,index=c.index)
    lowest_price_diet.append(result.fun)
d = {'Age/Sex Group': age_list, 'Lowest Price Diet ($ / Day)': lowest_price_diet}
total_dataframe = pd.DataFrame(data=d)
total_dataframe = total_dataframe.set_index("Age/Sex Group")
total_dataframe["Population Percentage"] = age_breakdown
total_dataframe["Weighted Cost"] = total_dataframe["Population Percentage"].astype(float) / 100 * total_dataframe["Lowest Price Diet ($ / Day)"]

print("The average weighted minimum cost diet is $" + str(np.round(sum(total_dataframe["Weighted Cost"]),2)) +" per person")
total_dataframe<jupyter_output>The average weighted minimum cost diet is $5.38 per person<jupyter_text>Regional Household Expenditures<jupyter_code>all_households = pd.read_csv('faps_household_puf.csv', dtype=str)
northeast = all_households[all_households["region"] == "1"]
midwest = all_households[all_households["region"] == "2"]
south = all_households[all_households["region"] == "3"]
west = all_households[all_households["region"] == "4"]
hh_choice = pd.read_csv('faps_fafhnutrient_puf.csv', dtype=str)
hh_choice["food"] = hh_choice["usdafoodcat4"].map(value_name_dict)
hh_choice['gramstotal'] = hh_choice['gramstotal'].astype(float)
northeast_hh = hh_choice[hh_choice["hhnum"].isin(northeast["hhnum"])]
midwest_hh = hh_choice[hh_choice["hhnum"].isin(midwest["hhnum"])]
south_hh = hh_choice[hh_choice["hhnum"].isin(south["hhnum"])]
west_hh = hh_choice[hh_choice["hhnum"].isin(west["hhnum"])]

northeast_hh['gramstotal'] = northeast_hh['gramstotal'].astype(float)
midwest_hh['gramstotal'] = midwest_hh['gramstotal'].astype(float)
south_hh['gramstotal'] = south_hh['gramstotal'].astype(float)
west_hh['gramstotal'] = west_hh['gramstotal'].astype(float)
import re
def numbers(value):
    return re.findall(r'\d+.\d+', value)
list_price = Prices.astype(str).apply(numbers)
def first_list(value):
    if len(value) >= 1:
        return float(value[0])
    else:
        return float(0)
corrected_prices = list_price.apply(first_list)
hh_stats = pd.read_csv('faps_individual_puf.csv', dtype=str)
hh_stats["hhnum"] = hh_stats["HHNUM"]
hh_stats_number_of_people = len(hh_stats)
hh_stats_average_number_days = np.average(hh_stats["INITFINALDAYS"].astype(float))
northeast_hh_number_of_people = len(hh_stats[hh_stats["HHNUM"].isin(northeast["hhnum"])])
midwest_hh_number_of_people = len(hh_stats[hh_stats["HHNUM"].isin(midwest["hhnum"])])
south_hh_number_of_people = len(hh_stats[hh_stats["HHNUM"].isin(south["hhnum"])])
west_hh_number_of_people = len(hh_stats[hh_stats["HHNUM"].isin(west["hhnum"])])
min_cost = ['filler']
regions = [hh_choice, northeast_hh, midwest_hh, south_hh, west_hh, min_cost]
values = []
households = []
number_of_people = []
monthly_income = []
measurement_time = []
for region in regions:
    if len(region) == 1:
        values.append(5.38)
        region = hh_choice
        food_df = pd.DataFrame()
        food_df["Foods"] = region.groupby("food")["gramstotal"].sum() / len(hh_stats[hh_stats["HHNUM"].isin(region["hhnum"])])
        monthly_income.append(np.average(all_households[all_households["hhnum"].isin(region["hhnum"])]['inchhavg_r'].astype(float)))
        households.append(len(hh_stats[hh_stats["hhnum"].isin(region["hhnum"])]["hhnum"].unique()))
        number_of_people.append(len(hh_stats[hh_stats["HHNUM"].isin(region["hhnum"])]))
        region['gramstotal'] = region['gramstotal'].astype(float)
        food_df["Price per Gram"] = corrected_prices / 100
        food_df["Amount Spent"] = food_df["Foods"] * food_df["Price per Gram"] 
        number_of_days = np.average(hh_stats[hh_stats["HHNUM"].isin(region["hhnum"])]["INITFINALDAYS"].astype(int))
        measurement_time.append(number_of_days)
        food_df["Amount Spent Daily"] = food_df["Amount Spent"] / number_of_days
    else:
        food_df = pd.DataFrame()
        food_df["Foods"] = region.groupby("food")["gramstotal"].sum() / len(hh_stats[hh_stats["HHNUM"].isin(region["hhnum"])])
        households.append(len(hh_stats[hh_stats["HHNUM"].isin(region["hhnum"])]["hhnum"].unique()))
        monthly_income.append(np.average(all_households[all_households["hhnum"].isin(region["hhnum"])]['inchhavg_r'].astype(float)))
        number_of_people.append(len(hh_stats[hh_stats["HHNUM"].isin(region["hhnum"])]))
        region['gramstotal'] = region['gramstotal'].astype(float)
        food_df["Price per Gram"] = corrected_prices / 100
        food_df["Amount Spent"] = food_df["Foods"] * food_df["Price per Gram"] 
        number_of_days = np.average(hh_stats[hh_stats["HHNUM"].isin(region["hhnum"])]["INITFINALDAYS"].astype(int))
        food_df["Amount Spent Daily"] = food_df["Amount Spent"] / number_of_days
        measurement_time.append(number_of_days)
        values.append(food_df["Amount Spent Daily"].sum())

values = np.round(values, 2)
regions = ["Total", "Northeast", "Midwest", "South", "West", "Total Min Cost"]
total_df = pd.DataFrame()
total_df["Region"] = regions
total_df["Amount Spent Daily ($)"] = values
total_df["Number of Households"] = households
total_df["Number of People"] = number_of_people
ratio = [int(p) / int(h) for p,h in zip(number_of_people, households)]
ratio = np.round(ratio, 2)
total_df["Average People per Household"] = ratio
#total_df["Average length of Interview"] = measurement_time
total_df["Daily Income ($)"] = [np.round(int(p) / 30, 2) for p in monthly_income]
total_df["Proportion of Income Spent on Food (%)"] =  total_df["Amount Spent Daily ($)"] / total_df["Daily Income ($)"] * 100
total_df= total_df.set_index('Region')
total_df

white = hh_stats[hh_stats["RACECAT_R"] == "1"]
black = hh_stats[hh_stats["RACECAT_R"] == "2"]
american_indian = hh_stats[hh_stats["RACECAT_R"] == "3"]
asian = hh_stats[hh_stats["RACECAT_R"] == "4"]
other = hh_stats[hh_stats["RACECAT_R"] == "5"]
multiple = hh_stats[hh_stats["RACECAT_R"] == "6"]

white_hh = hh_choice[hh_choice["hhnum"].isin(white["hhnum"])]
black_hh = hh_choice[hh_choice["hhnum"].isin(black["hhnum"])]
american_indian_hh = hh_choice[hh_choice["hhnum"].isin(american_indian["hhnum"])]
asian_hh = hh_choice[hh_choice["hhnum"].isin(asian["hhnum"])]
other_hh = hh_choice[hh_choice["hhnum"].isin(other["hhnum"])]
multiple_hh = hh_choice[hh_choice["hhnum"].isin(multiple["hhnum"])]

white_hh['gramstotal'] = white_hh['gramstotal'].astype(float)
black_hh['gramstotal'] = black_hh['gramstotal'].astype(float)
american_indian_hh['gramstotal'] = american_indian_hh['gramstotal'].astype(float)
asian_hh['gramstotal'] = asian_hh['gramstotal'].astype(float)
other_hh['gramstotal'] = other_hh['gramstotal'].astype(float)
multiple_hh['gramstotal'] = multiple_hh['gramstotal'].astype(float)
min_cost = ['filler']
races = [hh_choice, white_hh, black_hh, american_indian_hh, asian_hh, other_hh, multiple_hh]
values = []
households = []
monthly_income = []
number_of_people = []
measurement_time = []
for race in races:
    if len(region) == 1:
        values.append(1.38)
        race = hh_choice
        food_df = pd.DataFrame()
        food_df["Foods"] = race.groupby("food")["gramstotal"].sum() / len(hh_stats[hh_stats["HHNUM"].isin(race["hhnum"])])
        households.append(len(hh_stats[hh_stats["HHNUM"].isin(race["hhnum"])]["hhnum"].unique()))
        monthly_income.append(np.average(all_households[all_households["hhnum"].isin(race["hhnum"])]['inchhavg_r'].astype(float)))
        number_of_people.append(len(hh_stats[hh_stats["HHNUM"].isin(race["hhnum"])]))
        race['gramstotal'] = race['gramstotal'].astype(float)
        food_df["Price per Gram"] = corrected_prices / 100
        food_df["Amount Spent"] = food_df["Foods"] * food_df["Price per Gram"] 
        number_of_days = np.average(hh_stats[hh_stats["HHNUM"].isin(race["hhnum"])]["INITFINALDAYS"].astype(int))
        measurement_time.append(number_of_days)
        food_df["Amount Spent Daily"] = food_df["Amount Spent"] / number_of_days
    else:
        food_df = pd.DataFrame()
        food_df["Foods"] = race.groupby("food")["gramstotal"].sum() / len(hh_stats[hh_stats["HHNUM"].isin(race["hhnum"])])
        households.append(len(hh_stats[hh_stats["HHNUM"].isin(race["hhnum"])]["hhnum"].unique()))
        monthly_income.append(np.average(all_households[all_households["hhnum"].isin(race["hhnum"])]['inchhavg_r'].astype(float)))
        number_of_people.append(len(hh_stats[hh_stats["HHNUM"].isin(race["hhnum"])]))
        race['gramstotal'] = race['gramstotal'].astype(float)
        food_df["Price per Gram"] = corrected_prices / 100
        food_df["Amount Spent"] = food_df["Foods"] * food_df["Price per Gram"] 
        number_of_days = np.average(hh_stats[hh_stats["HHNUM"].isin(race["hhnum"])]["INITFINALDAYS"].astype(int))
        measurement_time.append(number_of_days)
        food_df["Amount Spent Daily"] = food_df["Amount Spent"] / number_of_days
        values.append(food_df["Amount Spent Daily"].sum())

values = np.round(values, 2)
races = ["Total", "White", "Black", "American Indian", "Asian", "Other", "Multiple"]
total_df = pd.DataFrame()
total_df["Race"] = races
total_df["Amount Spent Daily ($)"] = values
total_df["Number of Households"] = households
total_df["Number of People"] = number_of_people
ratio = [int(p) / int(h) for p,h in zip(number_of_people, households)]
ratio = np.round(ratio, 2)
total_df["Average People per Household"] = ratio
#total_df["Average length of Interview"] = measurement_time
total_df["Daily Income ($)"] = [np.round(int(p) / 30, 2) for p in monthly_income]
total_df["Proportion of Income Spent on Food (%)"] =  total_df["Amount Spent Daily ($)"] / total_df["Daily Income ($)"] * 100
total_df= total_df.set_index('Race')
total_df<jupyter_output><empty_output><jupyter_text>Nutritional Thresholds<jupyter_code>food_df = pd.DataFrame()
food_df["Amount in Hectograms"] = hh_choice.groupby("food")["gramstotal"].sum() / len(hh_stats) / 100
hh_choice['gramstotal'] = hh_choice['gramstotal'].astype(float)
food_df["Price per Hectogram"] = corrected_prices
food_df["Amount Spent"] = food_df["Amount in Hectograms"] * food_df["Price per Hectogram"] 
results = np.abs(b).align(D.fillna(0) @ food_df[["Amount in Hectograms"]]["Amount in Hectograms"], 'inner')[1]
recommendation = np.abs(b).align(D.fillna(0) @ food_df[["Amount in Hectograms"]]["Amount in Hectograms"], 'inner')[0]
tab = pd.DataFrame()
tab["Outcome"] = results
tab["Recommendations"] = recommendation
print("\nWith the following nutritional outcomes of interest (All Locations):")
print(tab)

fixed_D = D.loc[:,D.columns.to_series().isin(food_df.index.to_series()).tolist()]
food_df = pd.DataFrame()
food_df["Amount in Hectograms"] = northeast_hh.groupby("food")["gramstotal"].sum() / len(hh_stats[hh_stats["HHNUM"].isin(northeast_hh["hhnum"])]) / 100
northeast_hh['gramstotal'] = northeast_hh['gramstotal'].astype(float)
food_df["Price per Hectogram"] = corrected_prices
food_df["Amount Spent"] = food_df["Amount in Hectograms"] * food_df["Price per Hectogram"] 
fixed_D = D.loc[:,D.columns.to_series().isin(food_df.index.to_series()).tolist()]
results = np.abs(b).align(fixed_D.fillna(0) @ food_df[["Amount in Hectograms"]]["Amount in Hectograms"], 'inner')[1]
recommendation = np.abs(b).align(fixed_D.fillna(0) @ food_df[["Amount in Hectograms"]]["Amount in Hectograms"], 'inner')[0]
tab = pd.DataFrame()
tab["Outcome"] = results
tab["Recommendations"] = recommendation
print("\nWith the following nutritional outcomes of interest (Northeast):")
print(tab)
food_df = pd.DataFrame()
food_df["Amount in Hectograms"] = midwest_hh.groupby("food")["gramstotal"].sum() / len(hh_stats[hh_stats["HHNUM"].isin(midwest_hh["hhnum"])]) / 100
midwest_hh['gramstotal'] = midwest_hh['gramstotal'].astype(float)
food_df["Price per Hectogram"] = corrected_prices
food_df["Amount Spent"] = food_df["Amount in Hectograms"] * food_df["Price per Hectogram"] 
fixed_D = D.loc[:,D.columns.to_series().isin(food_df.index.to_series()).tolist()]
results = np.abs(b).align(fixed_D.fillna(0) @ food_df[["Amount in Hectograms"]]["Amount in Hectograms"], 'inner')[1]
recommendation = np.abs(b).align(fixed_D.fillna(0) @ food_df[["Amount in Hectograms"]]["Amount in Hectograms"], 'inner')[0]
tab = pd.DataFrame()
tab["Outcome"] = results
tab["Recommendations"] = recommendation
print("\nWith the following nutritional outcomes of interest (Midwest):")
print(tab)
food_df = pd.DataFrame()
food_df["Amount in Hectograms"] = south_hh.groupby("food")["gramstotal"].sum() / len(hh_stats[hh_stats["HHNUM"].isin(south_hh["hhnum"])]) / 100
south_hh['gramstotal'] = south_hh['gramstotal'].astype(float)
food_df["Price per Hectogram"] = corrected_prices
food_df["Amount Spent"] = food_df["Amount in Hectograms"] * food_df["Price per Hectogram"] 
fixed_D = D.loc[:,D.columns.to_series().isin(food_df.index.to_series()).tolist()]
results = np.abs(b).align(fixed_D.fillna(0) @ food_df[["Amount in Hectograms"]]["Amount in Hectograms"], 'inner')[1]
recommendation = np.abs(b).align(fixed_D.fillna(0) @ food_df[["Amount in Hectograms"]]["Amount in Hectograms"], 'inner')[0]
tab = pd.DataFrame()
tab["Outcome"] = results
tab["Recommendations"] = recommendation
print("\nWith the following nutritional outcomes of interest (South):")
print(tab)
food_df = pd.DataFrame()
food_df["Amount in Hectograms"] = west_hh.groupby("food")["gramstotal"].sum() / len(hh_stats[hh_stats["HHNUM"].isin(west_hh["hhnum"])]) / 100
west_hh['gramstotal'] = west_hh['gramstotal'].astype(float)
food_df["Price per Hectogram"] = corrected_prices
food_df["Amount Spent"] = food_df["Amount in Hectograms"] * food_df["Price per Hectogram"] 
fixed_D = D.loc[:,D.columns.to_series().isin(food_df.index.to_series()).tolist()]
results = np.abs(b).align(fixed_D.fillna(0) @ food_df[["Amount in Hectograms"]]["Amount in Hectograms"], 'inner')[1]
recommendation = np.abs(b).align(fixed_D.fillna(0) @ food_df[["Amount in Hectograms"]]["Amount in Hectograms"], 'inner')[0]
tab = pd.DataFrame()
tab["Outcome"] = results
tab["Recommendations"] = recommendation
print("\nWith the following nutritional outcomes of interest (West):")
print(tab)<jupyter_output>/srv/app/venv/lib/python3.6/site-packages/ipykernel_launcher.py:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  This is separate from the ipykernel package so we can avoid doing imports until<jupyter_text>Demand Estimation Cleaning<jupyter_code>hh_choice['gramstotal'] = hh_choice['gramstotal'].astype(float)
households = all_households["hhnum"].tolist()
regions = all_households["region"].tolist()
household_region_dict = dict(zip(households,regions))
hh_choice["region"] = hh_choice["hhnum"].map(household_region_dict)
grouped = hh_choice.groupby(["hhnum","region", "food"])["gramstotal"].sum().unstack('food')
grouped = grouped.reset_index().rename_axis(None, axis=1)
grouped.insert(1, "year", [2013] * len(grouped))
grouped = grouped.rename(index=str, columns={"hhnum": "j", "year": "t", "region": "m"})<jupyter_output><empty_output><jupyter_text>Household Data<jupyter_code>hh_stats = pd.read_csv('faps_individual_puf.csv', dtype=str)
hh_stats = pd.read_csv('faps_individual_puf.csv', dtype=str)
hh_stats = hh_stats[["HHNUM", "SEX", "AGE_R"]]
hh_stats = hh_stats[hh_stats["AGE_R"] != 'D']
hh_stats['AGE_R'] = hh_stats['AGE_R'].astype(float)
def men_function(x):
    if float(x["AGE_R"]) >= 18 and x["SEX"] == '1':
        return 1
    return 0
def boy_function(x):
    if float(x["AGE_R"]) < 18 and x["SEX"] == '1':
        return 1
    return 0
def women_function(x):
    if float(x["AGE_R"]) >= 18 and x["SEX"] == '2':
        return 1
    return 0
def girl_function(x):
    if float(x["AGE_R"]) < 18 and x["SEX"] == '2':
        return 1
    return 0
#hh_stats["ADULT"] = hh_stats["AGE_R"].apply(over_eighteen)
hh_stats["Men"] = hh_stats.apply(men_function, axis=1)
hh_stats["Women"] = hh_stats.apply(women_function, axis=1)
hh_stats["Boys"] = hh_stats.apply(boy_function, axis=1)
hh_stats["Girls"] = hh_stats.apply(girl_function, axis=1)


hh_stats.head(30)
def household_size(x):
    return np.log(x["Men"] + x["Women"] + x["Boys"] + x["Girls"])
hh_stats = hh_stats[["HHNUM", "Men", "Women", "Boys", "Girls"]]
hh_stats_grouped = hh_stats.groupby("HHNUM").sum().reset_index()
region_col = hh_stats_grouped["HHNUM"].map(household_region_dict)
hh_stats_grouped.insert(1, "region", region_col)
hh_stats_grouped.insert(1, "year", [2013] * len(hh_stats_grouped))
hh_stats_grouped["logHsize"] = hh_stats_grouped.apply(household_size, axis=1)
hh_stats_grouped = hh_stats_grouped.rename(index=str, columns={"HHNUM": "j", "year": "t", "region": "m"})
hh_stats_grouped
import pandas as pd
import numpy as np

x = grouped

# Change 'ICRISAT' to key of your own sheet in Sheets, above
z = hh_stats_grouped

# Assume a single market: (Comment this out to make each village its own market)
z['m'] = 1
x['m'] = 1

# Take logs of expenditures; call this y
y = np.log(x.replace(0,np.nan).set_index(['j','t','m']))

z.set_index(['j','t','m'],inplace=True)
z.head()
y.head()
import cfe

result = cfe.Result(y=y,z=z)

result<jupyter_output><empty_output><jupyter_text>First step of Estimation Recall that there are two steps to estimation; the first stepinvolves estimating the &ldquo;reduced form&rdquo; linear regression $$y_{it}^j = {a}_{it} + \delta_i'{z}^j_t + \epsilon_{it}^j.$$The Result class has code to estimate this in one line:<jupyter_code>result.get_reduced_form()<jupyter_output>/srv/app/venv/lib/python3.6/site-packages/pandas/core/frame.py:6211: FutureWarning: Sorting because non-concatenation axis is not aligned. A future version
of pandas will change to not sort by default.

To accept the future behavior, pass 'sort=False'.

To retain the current behavior and silence the warning, pass 'sort=True'.

  sort=sort)<jupyter_text>After running this we can examine the estimated coefficients $\delta$:<jupyter_code>result.delta.to_dataframe().unstack('k')<jupyter_output><empty_output><jupyter_text>Also the good-time constants $a_{it}$:<jupyter_code>result.a.to_dataframe().unstack('i')<jupyter_output><empty_output><jupyter_text>Second step of Estimation The second step involves using Singular Value Decomposition to findthe rank one matrix that best approximates the residuals $e_{it}^j$.This can be interpreted as$$    -\beta_i\log\lambda^j_t,$$where the $\log\lambda^j_t$ is the log of the marginal utility ofexpenditures (MUE) for household $j$ at time $t$, and where $\beta_i$ arethe corresponding &ldquo;Frisch elasticities&rdquo; that tell us how muchdemand changes as the MUE falls.Estimates can also be computed as a one-liner:<jupyter_code>result.get_beta().to_dataframe().sort_values(by=["beta"], ascending=False).head(10)<jupyter_output><empty_output><jupyter_text>That&rsquo;s all there is to estimation!  Note that we didn&rsquo;t estimatedemands for all goods&x2014;lots of goods didn&rsquo;t have enough observations,and were automatically dropped.  (This can be controlled using the`min_proportion_items` and `min_xproducts` attributes when oneinstantiates the result object.)Now, let&rsquo;s see how we did, by comparing total expenditures predicted by themodel we&rsquo;ve estimated with actual total expenditures:<jupyter_code>import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib.cm as cm

result.get_predicted_log_expenditures()

xbar = np.exp(result.y).sum(['m','i']).to_dataframe('xbar').replace(0,np.nan).squeeze()

xhat = result.get_predicted_expenditures().sum(['m','i']).to_dataframe('xhat').replace(0,np.nan).squeeze()

# Make dataframe of actual & predicted
df = pd.DataFrame({'Actual':np.log(xbar),'Predicted':np.log(xhat)})

df.plot.scatter(x='Predicted',y='Actual')<jupyter_output><empty_output><jupyter_text>Holding Households Constant<jupyter_code>house_result = result
house_result['z'] = house_result['z'] * 0
house_result.get_predicted_log_expenditures()

x = house_result.y.sum('i').to_dataframe().squeeze()
xhat = house_result.yhat.sum('i').to_dataframe().squeeze()

xbar = np.exp(result.y).sum(['m','i']).to_dataframe('xbar').replace(0,np.nan).squeeze()

xhat = house_result.get_predicted_expenditures().sum(['m','i']).to_dataframe('xhat').replace(0,np.nan).squeeze()

# Make dataframe of actual & predicted
df = pd.DataFrame({'Actual':np.log(xbar),'Predicted':np.log(xhat)})

df.plot.scatter(x='Predicted',y='Actual')<jupyter_output><empty_output><jupyter_text>Engle's Law<jupyter_code>result = cfe.Result(y=y,z=z)
result.get_predicted_log_expenditures()

cheese=result.get_predicted_expenditures().sel(i='Cheese',t=2013,m=1)
xbar = result.get_predicted_expenditures().sel(t=2013,m=1).sum('i')

wcheese = cheese/xbar

# Make dataframe of actual & predicted
df = pd.DataFrame({'Actual Expenditures':np.log(xbar),'Predicted Consumption of Cheese':wcheese})

df.plot.scatter(x='Actual Expenditures',y='Predicted Consumption of Cheese')
result.get_predicted_log_expenditures()

bread=result.get_predicted_expenditures().sel(i='Yeast breads',t=2013,m=1)
xbar = result.get_predicted_expenditures().sel(t=2013,m=1).sum('i')

wbread = bread/xbar

# Make dataframe of actual & predicted
df = pd.DataFrame({'Actual Expenditures':np.log(xbar),'Predicted Consumption of Bread':wbread})

df.plot.scatter(x='Actual Expenditures',y='Predicted Consumption of Bread')
result.get_predicted_log_expenditures()

rice=result.get_predicted_expenditures().sel(i='Rice',t=2013,m=1)
xbar = result.get_predicted_expenditures().sel(t=2013,m=1).sum('i')

wrice = rice/xbar

# Make dataframe of actual & predicted
df = pd.DataFrame({'Actual Expenditures':np.log(xbar),'Predicted Consumption of Rice':wrice})

df.plot.scatter(x='Actual Expenditures',y='Predicted Consumption of Rice')
result.get_predicted_log_expenditures()

beef=result.get_predicted_expenditures().sel(i='Beef, excludes ground',t=2013,m=1)
xbar = result.get_predicted_expenditures().sel(t=2013,m=1).sum('i')

wbeef = beef/xbar

# Make dataframe of actual & predicted
df = pd.DataFrame({'Actual Expenditures':np.log(xbar),'Predicted Consumption of Beef, excludes Ground':wbeef})

df.plot.scatter(x='Actual Expenditures',y='Predicted Consumption of Beef, excludes Ground')<jupyter_output>/srv/app/venv/lib/python3.6/site-packages/xarray/core/computation.py:561: RuntimeWarning: divide by zero encountered in log
  result_data = func(*input_data)
{'fraction_non_alphanumeric': 0.17010441472797216, 'fraction_numerical': 0.012676314343286315, 'mean_word_length': 6.06235446313066, 'pattern_counts': {'":': 8, '<': 50, '<?xml version=': 0, '>': 50, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '114215', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 10380, 'n_tokens_neox': 9830, 'n_words': 1839, 'path': 'Project 3.ipynb', 'repo_name': 'nathanorwig/Atwater-Min-Cost-Diet', 'size': '27330'}