import numpy as np
import datetime
import pandas as pd
import glob
import matplotlib.pylab as plt
import sys
import statsmodels.api as sm
import warnings
warnings.filterwarnings("ignore")
import os
import csv

# Weighted percentile calculator
def wquantile(x,q):           
	xsort = x.sort_values(x.columns[0])
	xsort['index'] = range(len(x))
	p = q * x[x.columns[1]].sum()
	pop = float(xsort[xsort.columns[1]][xsort['index']==0])
	i = 0
	while pop < p:
		pop = pop + float(xsort[xsort.columns[1]][xsort['index']==i+1])
		i = i + 1
	return xsort[xsort.columns[0]][xsort['index']==i]



################################################################
################################################################
####
#### This program file is part of the article "What Do Survey Data Tell Us about U.S. Businesses?"
#### It is used for generating moments from SIPP for each SIPP panel. 
####  	
################################################################
################################################################

#sipp_panel     = int(input("Please enter SIPP Panel Year (1996, 2001, 2004, or 2008): "))

for sipp_panel in [1996, 2001, 2004, 2008]:

	print("--- GENERATING MOMENTS FOR SIPP PANEL %s ---" % sipp_panel)

	sipp_panel_str = str(sipp_panel)
	sipp_merge     = pd.read_csv('sipp_'+sipp_panel_str+'.csv')

	max_year 		= np.max(sipp_merge.year)
	min_year 		= np.min(sipp_merge.year)
	min_month = 1
	max_month = 12

	# Create csv for results if it doesn't exist
	sipp_dir = os.getcwd()+'/' # Be sure it's the same directory where the sipp_xxxx.csv files are.
	if os.path.isfile(sipp_dir+'ta11.csv')==0:
		my_dataframe = pd.DataFrame()
		col_names = ["Tax Year", "SP VW_Mean", "SP EW_Mean", "SP p25", "SP p50", "SP p75", "P VW_Mean", "P EW_Mean", "P p25", "P p50", "P p75", "U VW_Mean", "U EW_Mean", "U p25", "U p50", "U p75"]
		my_dataframe = pd.DataFrame(columns=col_names)
		my_dataframe = my_dataframe.set_index('Tax Year')
		my_dataframe.to_csv(sipp_dir+'ta11.csv', index=True)

	###### IRS ######

	# Partnerships

	# 1) Number of partnerships over time - Topical, so sum across months

	# Sum ownership shares of those who report having a partnership (weighted).

	partnership_num = np.zeros((max_year - min_year  + 1, max_month-min_month+1))
	partnership_net_income_num = np.zeros((max_year - min_year  + 1, max_month-min_month+1))
	partnership_net_loss_num = np.zeros((max_year - min_year  + 1, max_month-min_month+1))

	partnership_zero_ownership = np.zeros((max_year - min_year  + 1, max_month-min_month+1))

	partnership_select_1 	= sipp_merge['LEG_SOLE_PART_1'] == 2
	partnership_select_2 	= sipp_merge['LEG_SOLE_PART_2'] == 2

	gain_select_1 			= sipp_merge['NET_PFT_LOSS_BUS_1'] > 0
	gain_select_2 			= sipp_merge['NET_PFT_LOSS_BUS_2'] > 0

	loss_select_1 			= sipp_merge['NET_PFT_LOSS_BUS_1'] <= 0
	loss_select_2 			= sipp_merge['NET_PFT_LOSS_BUS_2'] <= 0

	for iy in np.arange(min_year,max_year+1):
		for im in np.arange(min_month,max_month+1):
			year_select = sipp_merge['year'] == iy
			month_select = sipp_merge['month'] == im

			partnership_num_1 = np.sum(sipp_merge[year_select & month_select & partnership_select_1].PCT_OWN_BUS_1/100 * sipp_merge[month_select & year_select & partnership_select_1].PER_WEIGHT)
			partnership_num_2 = np.sum(sipp_merge[year_select & month_select & partnership_select_2].PCT_OWN_BUS_2/100 * sipp_merge[month_select & year_select & partnership_select_2].PER_WEIGHT)

			partnership_net_income_num_1 = np.sum(sipp_merge[year_select & month_select & partnership_select_1 & gain_select_1].PCT_OWN_BUS_1/100 * sipp_merge[month_select & year_select & partnership_select_1 & gain_select_1].PER_WEIGHT)
			partnership_net_income_num_2 = np.sum(sipp_merge[year_select & month_select & partnership_select_2 & gain_select_2].PCT_OWN_BUS_2/100 * sipp_merge[month_select & year_select & partnership_select_2 & gain_select_2].PER_WEIGHT)

			partnership_net_loss_num_1 = np.sum(sipp_merge[year_select & month_select & partnership_select_1 & loss_select_1].PCT_OWN_BUS_1/100 * sipp_merge[month_select & year_select & partnership_select_1 & loss_select_1].PER_WEIGHT)
			partnership_net_loss_num_2 = np.sum(sipp_merge[year_select & month_select & partnership_select_2 & loss_select_2].PCT_OWN_BUS_2/100 * sipp_merge[month_select & year_select & partnership_select_2 & loss_select_2].PER_WEIGHT)

			partnership_num[iy-min_year,im-1] 				= partnership_num_1 + partnership_num_2
			partnership_net_income_num[iy-min_year,im-1] 	= partnership_net_income_num_1 + partnership_net_income_num_2
			partnership_net_loss_num[iy-min_year,im-1] 		= partnership_net_loss_num_1 + partnership_net_loss_num_2


	if sipp_panel == 1996: # Topical module time change
		partnership_num_annual 				= np.zeros(5)
		partnership_num_annual[1] 			= np.sum(partnership_num[1,10:])/10**3 + np.sum(partnership_num[2,0:2])/10**3
		partnership_num_annual[2] 			= np.sum(partnership_num[2,10:])/10**3 + np.sum(partnership_num[3,0:2])/10**3
		partnership_num_annual[3] 			= np.sum(partnership_num[3,8:])/10**3
		partnership_num_annual[4] 			= np.sum(partnership_num[4,8:])/10**3

		partnership_net_income_num_annual 	 = np.zeros(5)
		partnership_net_income_num_annual[1] = np.sum(partnership_net_income_num[1,10:])/10**3 + np.sum(partnership_net_income_num[2,0:2])/10**3
		partnership_net_income_num_annual[2] = np.sum(partnership_net_income_num[2,10:])/10**3 + np.sum(partnership_net_income_num[3,0:2])/10**3
		partnership_net_income_num_annual[3] = np.sum(partnership_net_income_num[3,8:])/10**3
		partnership_net_income_num_annual[4] = np.sum(partnership_net_income_num[4,8:])/10**3
		
		partnership_net_loss_num_annual 	= np.zeros(5)
		partnership_net_loss_num_annual[1] 	= np.sum(partnership_net_loss_num[1,10:])/10**3 + np.sum(partnership_net_loss_num[2,0:2])/10**3
		partnership_net_loss_num_annual[2] 	= np.sum(partnership_net_loss_num[2,10:])/10**3 + np.sum(partnership_net_loss_num[3,0:2])/10**3
		partnership_net_loss_num_annual[3] 	= np.sum(partnership_net_loss_num[3,8:])/10**3
		partnership_net_loss_num_annual[4] 	= np.sum(partnership_net_loss_num[4,8:])/10**3
	else:
		partnership_num_annual 				= np.sum(partnership_num, axis=1)/10**3
		partnership_net_income_num_annual 	= np.sum(partnership_net_income_num, axis=1)/10**3
		partnership_net_loss_num_annual 	= np.sum(partnership_net_loss_num, axis=1)/10**3

	# 2) Number of partners over time - Core, counting number, so we should average across months with complete info
	# Sum those who report owning at least one partnership.

	partner_num = np.zeros((max_year - min_year  + 1, max_month-min_month+1))

	sipp_merge['partnership_own_1_temp'] 	= sipp_merge['LEG_SOLE_PART_1'] == 2
	sipp_merge['partnership_own_2_temp'] 	= sipp_merge['LEG_SOLE_PART_2'] == 2
	#sipp_merge['partnership_own_temp'] 		= sipp_merge[['partnership_own_1_temp','partnership_own_2_temp']].max(axis=1) # 2 partnerships, coutns as 1 partner
	sipp_merge['partnership_own_temp'] 		= sipp_merge[['partnership_own_1_temp','partnership_own_2_temp']].sum(axis=1) # 2 partnerships, coutns as 2 partners

	for iy in np.arange(min_year,max_year+1):
		for im in np.arange(min_month,max_month+1):
			year_select = sipp_merge['year'] == iy
			month_select = sipp_merge['month'] == im

			partner_num[iy-min_year,im-1] = np.sum(sipp_merge[year_select & month_select].partnership_own_temp * sipp_merge[month_select & year_select].PER_WEIGHT)

	indicator_mat = np.ones((max_year - min_year  + 1, max_month-min_month+1))
	min_year_drop = np.sum(partner_num[0,:] == 0) + 3
	max_year_drop = np.minimum(np.sum(partner_num[-1,:] == 0) + 3, 12)
	indicator_mat[0,0:min_year_drop] = 0
	indicator_mat[-1,max_month-max_year_drop:] = 0

	partner_num_annual = np.sum(partner_num*indicator_mat, axis=1)/np.sum(indicator_mat, axis=1)/10**6

	sipp_merge.drop(['partnership_own_1_temp', 'partnership_own_2_temp', 'partnership_own_temp'], axis=1)


	if sipp_panel >= 2004:
		partnership_net_income = np.zeros((max_year - min_year  + 1, max_month-min_month+1))
		partnership_net_income_amt = np.zeros((max_year - min_year  + 1, max_month-min_month+1))
		partnership_net_loss_amt = np.zeros((max_year - min_year  + 1, max_month-min_month+1))

		partnership_select_1 = sipp_merge['LEG_SOLE_PART_1'] == 2
		partnership_select_2 = sipp_merge['LEG_SOLE_PART_2'] == 2

		gain_select_1 		 	= sipp_merge['NET_PFT_LOSS_BUS_1'] > 0
		gain_select_2 			= sipp_merge['NET_PFT_LOSS_BUS_2'] > 0

		loss_select_1 		 	= sipp_merge['NET_PFT_LOSS_BUS_1'] <= 0
		loss_select_2 			= sipp_merge['NET_PFT_LOSS_BUS_2'] <= 0

		for iy in np.arange(min_year,max_year+1):
			for im in np.arange(min_month,max_month+1):

				year_select = sipp_merge['year'] == iy
				month_select = sipp_merge['month'] == im

				partnership_net_income_1 		= np.sum(sipp_merge[year_select & month_select & partnership_select_1].NET_PFT_LOSS_BUS_1 * sipp_merge[month_select & year_select & partnership_select_1].PER_WEIGHT)
				partnership_net_income_2 		= np.sum(sipp_merge[year_select & month_select & partnership_select_2].NET_PFT_LOSS_BUS_2 * sipp_merge[month_select & year_select & partnership_select_2].PER_WEIGHT)

				partnership_net_income_amt_1 	= np.sum(sipp_merge[year_select & month_select & partnership_select_1 & gain_select_1].NET_PFT_LOSS_BUS_1 * sipp_merge[month_select & year_select & partnership_select_1 & gain_select_1].PER_WEIGHT)
				partnership_net_income_amt_2 	= np.sum(sipp_merge[year_select & month_select & partnership_select_2 & gain_select_2].NET_PFT_LOSS_BUS_2 * sipp_merge[month_select & year_select & partnership_select_2 & gain_select_2].PER_WEIGHT)

				partnership_net_loss_amt_1 		= np.sum(sipp_merge[year_select & month_select & partnership_select_1 & loss_select_1].NET_PFT_LOSS_BUS_1 * sipp_merge[month_select & year_select & partnership_select_1 & loss_select_1].PER_WEIGHT)
				partnership_net_loss_amt_2 		= np.sum(sipp_merge[year_select & month_select & partnership_select_2 & loss_select_2].NET_PFT_LOSS_BUS_2 * sipp_merge[month_select & year_select & partnership_select_2 & loss_select_2].PER_WEIGHT)

				partnership_net_income[iy-min_year,im-1] 		= partnership_net_income_1 + partnership_net_income_2 			# sum over each year with complete info
				partnership_net_income_amt[iy-min_year,im-1] 	= partnership_net_income_amt_1 + partnership_net_income_amt_2 	# sum over each year with complete info
				partnership_net_loss_amt[iy-min_year,im-1] 		= partnership_net_loss_amt_1 + partnership_net_loss_amt_2 		# sum over each year with complete info


		partnership_net_income_annual 			= np.sum(partnership_net_income, axis=1)/10**9 # sum income across 12 months, divide by 4 bec value is for ref period
		partnership_net_income_annual[0] 		= 0
		partnership_net_income_annual[-1] 		= 0
		partnership_net_income_annual 			= partnership_net_income_annual/4

		partnership_net_income_amt_annual 		= np.sum(partnership_net_income_amt, axis=1)/10**9 # sum income across 12 months, divide by 4 bec value is for ref period
		partnership_net_income_amt_annual[0]  	= 0
		partnership_net_income_amt_annual[-1] 	= 0
		partnership_net_income_amt_annual 		= partnership_net_income_amt_annual/4

		partnership_net_loss_amt_annual 		= np.sum(partnership_net_loss_amt, axis=1)/10**9 # sum income across 12 months, divide by 4 bec value is for ref period
		partnership_net_loss_amt_annual[0]  	= 0
		partnership_net_loss_amt_annual[-1] 	= 0
		partnership_net_loss_amt_annual 		= partnership_net_loss_amt_annual/4

	elif sipp_panel < 2004:

		#####
		sipp_merge['PCT_OWN_BUS_1'].loc[np.isnan(sipp_merge['PCT_OWN_BUS_1'])] = 0
		sipp_merge['PCT_OWN_BUS_2'].loc[np.isnan(sipp_merge['PCT_OWN_BUS_2'])] = 0
		sipp_merge['PCT_OWN_BUS_1_COPY'] = 0
		sipp_merge['PCT_OWN_BUS_2_COPY'] = 0

		grouped = sipp_merge.groupby(['PID','year']) # similar to bysort in STATA

		max_share = lambda x: x.max()
		sipp_merge['PCT_OWN_BUS_1_COPY'] = (sipp_merge.groupby(['PID','year'])['PCT_OWN_BUS_1'].transform(max_share))
		sipp_merge['PCT_OWN_BUS_2_COPY'] = (sipp_merge.groupby(['PID','year'])['PCT_OWN_BUS_2'].transform(max_share))

		#####
		partnership_net_income = np.zeros((max_year - min_year  + 1, max_month-min_month+1))
		partnership_net_income_amt = np.zeros((max_year - min_year  + 1, max_month-min_month+1))
		partnership_net_loss_amt = np.zeros((max_year - min_year  + 1, max_month-min_month+1))

		partnership_select_1 = sipp_merge['LEG_SOLE_PART_1'] == 2
		partnership_select_2 = sipp_merge['LEG_SOLE_PART_2'] == 2

		gain_select_1 		 	= sipp_merge['NET_PFT_LOSS_BUS_1'] > 0
		gain_select_2 			= sipp_merge['NET_PFT_LOSS_BUS_2'] > 0

		loss_select_1 		 	= sipp_merge['NET_PFT_LOSS_BUS_1'] <= 0
		loss_select_2 			= sipp_merge['NET_PFT_LOSS_BUS_2'] <= 0

		for iy in np.arange(min_year,max_year+1):
			for im in np.arange(min_month,max_month+1):

				year_select = sipp_merge['year'] == iy
				month_select = sipp_merge['month'] == im

				partnership_net_income_1 		= np.sum(sipp_merge[year_select & month_select & partnership_select_1].NET_PFT_LOSS_BUS_1 * sipp_merge[month_select & year_select & partnership_select_1].PER_WEIGHT)
				partnership_net_income_2 		= np.sum(sipp_merge[year_select & month_select & partnership_select_2].NET_PFT_LOSS_BUS_2 * sipp_merge[month_select & year_select & partnership_select_2].PER_WEIGHT)

				partnership_net_income_amt_1 	= np.sum(sipp_merge[year_select & month_select & partnership_select_1 & gain_select_1].NET_PFT_LOSS_BUS_1 * sipp_merge[month_select & year_select & partnership_select_1 & gain_select_1].PER_WEIGHT)
				partnership_net_income_amt_2 	= np.sum(sipp_merge[year_select & month_select & partnership_select_2 & gain_select_2].NET_PFT_LOSS_BUS_2 * sipp_merge[month_select & year_select & partnership_select_2 & gain_select_2].PER_WEIGHT)

				partnership_net_loss_amt_1 		= np.sum(sipp_merge[year_select & month_select & partnership_select_1 & loss_select_1].NET_PFT_LOSS_BUS_1 * sipp_merge[month_select & year_select & partnership_select_1 & loss_select_1].PER_WEIGHT)
				partnership_net_loss_amt_2 		= np.sum(sipp_merge[year_select & month_select & partnership_select_2 & loss_select_2].NET_PFT_LOSS_BUS_2 * sipp_merge[month_select & year_select & partnership_select_2 & loss_select_2].PER_WEIGHT)

				partnership_net_income[iy-min_year,im-1] 		= partnership_net_income_1 + partnership_net_income_2 			# sum over each year with complete info
				partnership_net_income_amt[iy-min_year,im-1] 	= partnership_net_income_amt_1 + partnership_net_income_amt_2 	# sum over each year with complete info
				partnership_net_loss_amt[iy-min_year,im-1] 		= partnership_net_loss_amt_1 + partnership_net_loss_amt_2 		# sum over each year with complete info


		partnership_net_income_annual 			= np.sum(partnership_net_income, axis=1)/10**9 # sum income across 12 months, divide by 4 bec value is for ref period
		partnership_net_income_annual[0] 		= 0
		partnership_net_income_annual[1] 		= 0
		partnership_net_income_annual[-1] 		= 0
		partnership_net_income_annual 			= partnership_net_income_annual/4

		partnership_net_income_amt_annual 		= np.sum(partnership_net_income_amt, axis=1)/10**9 # sum income across 12 months, divide by 4 bec value is for ref period
		partnership_net_income_amt_annual[0]  	= 0
		partnership_net_income_annual[1] 		= 0
		partnership_net_income_amt_annual[-1] 	= 0
		partnership_net_income_amt_annual 		= partnership_net_income_amt_annual/4

		partnership_net_loss_amt_annual 		= np.sum(partnership_net_loss_amt, axis=1)/10**9 # sum income across 12 months, divide by 4 bec value is for ref period
		partnership_net_loss_amt_annual[0]  	= 0
		partnership_net_income_annual[1] 		= 0
		partnership_net_loss_amt_annual[-1] 	= 0
		partnership_net_loss_amt_annual 		= partnership_net_loss_amt_annual/4




	# Sole Proprietorships

	# 1a) Number of returns, all, with net income, with net loss - Core, counting number, so we should average across months with complete info (corected income - take ownership as truth)
	# Sum those who report having at least one sole-prop. 
	# ASSUMPTION: We are assuming that those who have multiple sole-props file one income tax return with many Schedule C's to separately report business activity.

	sole_proprietor_num = np.zeros((max_year - min_year  + 1, max_month-min_month+1))
	sole_proprietor_net_income_num 	= np.zeros((max_year - min_year  + 1, max_month-min_month+1))
	sole_proprietor_net_loss_num 	= np.zeros((max_year - min_year  + 1, max_month-min_month+1))

	# Determine ownership of non-farm sole proprietorship (note that businessess can only be non-farm in SIPP, as inferred from the choices of industries)
	sipp_merge['sole_proprietorship_own_1_temp'] 	= (sipp_merge['LEG_SOLE_PART_1'] == 1)
	sipp_merge['sole_proprietorship_own_2_temp'] 	= (sipp_merge['LEG_SOLE_PART_2'] == 1)
	sipp_merge['sole_proprietorship_own_temp'] 		= sipp_merge[['sole_proprietorship_own_1_temp','sole_proprietorship_own_2_temp']].max(axis=1)

	# Calculate total net income/loss from sole-proprietorships, sum if respondent had two sole-props
	sipp_merge['sole_prop_net_pft_loss_temp'] = 0 
	sipp_merge['sole_prop_net_pft_loss_temp'].loc[sipp_merge.LEG_SOLE_PART_1 == 1] = sipp_merge.NET_PFT_LOSS_BUS_1
	sipp_merge['sole_prop_net_pft_loss_temp'].loc[sipp_merge.LEG_SOLE_PART_2 == 1] = sipp_merge.NET_PFT_LOSS_BUS_2
	sipp_merge['sole_prop_net_pft_loss_temp'].loc[(sipp_merge.LEG_SOLE_PART_1 == 1) & (sipp_merge.LEG_SOLE_PART_2 == 1)] = sipp_merge.NET_PFT_LOSS_BUS_1 + sipp_merge.NET_PFT_LOSS_BUS_2

	# Select gainers/losers
	gain_select = sipp_merge['sole_prop_net_pft_loss_temp'] > 0
	loss_select = sipp_merge['sole_prop_net_pft_loss_temp'] <= 0

	for iy in np.arange(min_year,max_year+1):
		for im in np.arange(min_month,max_month+1):
			year_select = sipp_merge['year'] == iy
			month_select = sipp_merge['month'] == im

			sole_proprietor_num 			[iy-min_year,im-1] 	= np.sum(sipp_merge[year_select & month_select].sole_proprietorship_own_temp * sipp_merge[month_select & year_select].PER_WEIGHT)
			sole_proprietor_net_income_num 	[iy-min_year,im-1] 	= np.sum(sipp_merge[year_select & month_select & gain_select].sole_proprietorship_own_temp * sipp_merge[month_select & year_select & gain_select].PER_WEIGHT)
			sole_proprietor_net_loss_num	[iy-min_year,im-1] 	= np.sum(sipp_merge[year_select & month_select & loss_select].sole_proprietorship_own_temp * sipp_merge[month_select & year_select & loss_select].PER_WEIGHT)

	indicator_mat = np.ones((max_year - min_year  + 1, max_month-min_month+1))
	min_year_drop = np.sum(sole_proprietor_num[0,:] == 0) + 3
	max_year_drop = np.minimum(np.sum(sole_proprietor_num[-1,:] == 0) + 3, 12)
	indicator_mat[0,0:min_year_drop] = 0
	indicator_mat[-1,max_month-max_year_drop:] = 0

	sole_proprietor_num_annual 				= np.sum(sole_proprietor_num*indicator_mat, axis=1)/np.sum(indicator_mat, axis=1)/10**3
	sole_proprietor_net_income_num_annual 	= np.sum(sole_proprietor_net_income_num*indicator_mat, axis=1)/np.sum(indicator_mat, axis=1)/10**3
	sole_proprietor_net_loss_num_annual 	= np.sum(sole_proprietor_net_loss_num*indicator_mat, axis=1)/np.sum(indicator_mat, axis=1)/10**3

	#sipp_merge.drop(['sole_proprietorship_own_1_temp', 'sole_proprietorship_own_2_temp', 'sole_proprietorship_own_temp', 'sole_prop_net_pft_loss_temp'], axis=1)


	# 2) Total net income or loss over time - Core, flow monetary, so sum across years with complete info

	# OPTION 1 - Use business's income x share

	if sipp_panel >= 2004:
		sole_proprietorship_net_income 		= np.zeros((max_year - min_year  + 1, max_month-min_month+1))
		sole_proprietorship_net_income_amt 	= np.zeros((max_year - min_year  + 1, max_month-min_month+1))
		sole_proprietorship_net_loss_amt 	= np.zeros((max_year - min_year  + 1, max_month-min_month+1))

		sole_proprietorship_select_1 = sipp_merge['LEG_SOLE_PART_1'] == 1
		sole_proprietorship_select_2 = sipp_merge['LEG_SOLE_PART_2'] == 1

		gain_select_1 		 	= sipp_merge['NET_PFT_LOSS_BUS_1'] > 0
		gain_select_2 			= sipp_merge['NET_PFT_LOSS_BUS_2'] > 0
		loss_select_1 		 	= sipp_merge['NET_PFT_LOSS_BUS_1'] <= 0
		loss_select_2 			= sipp_merge['NET_PFT_LOSS_BUS_2'] <= 0

		for iy in np.arange(min_year,max_year+1):
			for im in np.arange(min_month,max_month+1):

				year_select = sipp_merge['year'] == iy
				month_select = sipp_merge['month'] == im

				sole_proprietorship_net_income_1 		= np.sum(sipp_merge[year_select & month_select & sole_proprietorship_select_1].NET_PFT_LOSS_BUS_1 * sipp_merge[month_select & year_select & sole_proprietorship_select_1].PER_WEIGHT)
				sole_proprietorship_net_income_2 		= np.sum(sipp_merge[year_select & month_select & sole_proprietorship_select_2].NET_PFT_LOSS_BUS_2 * sipp_merge[month_select & year_select & sole_proprietorship_select_2].PER_WEIGHT)

				sole_proprietorship_net_income_amt_1 	= np.sum(sipp_merge[year_select & month_select & sole_proprietorship_select_1 & gain_select_1].NET_PFT_LOSS_BUS_1 * sipp_merge[month_select & year_select & sole_proprietorship_select_1 & gain_select_1].PER_WEIGHT)
				sole_proprietorship_net_income_amt_2 	= np.sum(sipp_merge[year_select & month_select & sole_proprietorship_select_2 & gain_select_2].NET_PFT_LOSS_BUS_2 * sipp_merge[month_select & year_select & sole_proprietorship_select_2 & gain_select_2].PER_WEIGHT)

				sole_proprietorship_net_loss_amt_1 		= np.sum(sipp_merge[year_select & month_select & sole_proprietorship_select_1 & loss_select_1].NET_PFT_LOSS_BUS_1 * sipp_merge[month_select & year_select & sole_proprietorship_select_1 & loss_select_1].PER_WEIGHT)
				sole_proprietorship_net_loss_amt_2 		= np.sum(sipp_merge[year_select & month_select & sole_proprietorship_select_2 & loss_select_2].NET_PFT_LOSS_BUS_2 * sipp_merge[month_select & year_select & sole_proprietorship_select_2 & loss_select_2].PER_WEIGHT)


				sole_proprietorship_net_income 		[iy-min_year,im-1] = sole_proprietorship_net_income_1 + sole_proprietorship_net_income_2 # sum over each year with complete info
				sole_proprietorship_net_income_amt 	[iy-min_year,im-1] = sole_proprietorship_net_income_amt_1 + sole_proprietorship_net_income_amt_2
				sole_proprietorship_net_loss_amt 	[iy-min_year,im-1] = sole_proprietorship_net_loss_amt_1 + sole_proprietorship_net_loss_amt_2

		sole_proprietorship_net_income_annual = np.sum(sole_proprietorship_net_income, axis=1)/10**9
		sole_proprietorship_net_income_annual[0] 	= 0
		sole_proprietorship_net_income_annual[-1] 	= 0
		sole_proprietorship_net_income_annual = sole_proprietorship_net_income_annual/4

		sole_proprietorship_net_income_amt_annual = np.sum(sole_proprietorship_net_income_amt, axis=1)/10**9
		sole_proprietorship_net_income_amt_annual[0] 	= 0
		sole_proprietorship_net_income_amt_annual[-1] 	= 0
		sole_proprietorship_net_income_amt_annual = sole_proprietorship_net_income_amt_annual/4

		sole_proprietorship_net_loss_amt_annual = np.sum(sole_proprietorship_net_loss_amt, axis=1)/10**9
		sole_proprietorship_net_loss_amt_annual[0] 	= 0
		sole_proprietorship_net_loss_amt_annual[-1] = 0
		sole_proprietorship_net_loss_amt_annual = sole_proprietorship_net_loss_amt_annual/4

	#sipp_merge.drop(['sole_proprietorship_net_income_1_temp', 'sole_proprietorship_net_income_2_temp'], axis=1)



	# 2) Sole prop and partnership income over time - Core, flow monetary, so sum across years with complete info 
	# SIPP, non-farm, unadjusted/adjusted

	if sipp_panel >= 2004:
		# SAMPLE 
		sipp_merge['PROP_INC_1'] = sipp_merge.INC_BUS_1*(sipp_merge.LEG_SOLE_PART_1>0)
		sipp_merge['PROP_INC_2'] = sipp_merge.INC_BUS_2*(sipp_merge.LEG_SOLE_PART_2>0)

		month_select 	= sipp_merge['month'] == 1
		year_select 	= sipp_merge['year'] == 2009
		prop_select_1 	= sipp_merge['LEG_SOLE_PART_1'] > 0
		prop_select_2 	= sipp_merge['LEG_SOLE_PART_2'] > 0

		# LOOP
		# Note that we do not weight because the question asked is how much income was derived from business, not what is total business income?

		prop_inc 		= np.zeros((max_year - min_year  + 1, max_month-min_month+1))

		prop_select_1 = sipp_merge['LEG_SOLE_PART_1'] > 0
		prop_select_2 = sipp_merge['LEG_SOLE_PART_2'] > 0

		for iy in np.arange(min_year,max_year+1):
			for im in np.arange(min_month,max_month+1):
				year_select = sipp_merge['year'] == iy
				month_select = sipp_merge['month'] == im

				# OPTION 1 
				#prop_inc_1 = np.sum(sipp_merge[year_select & month_select & prop_select_1].PROP_INC_1 * sipp_merge[month_select & year_select & prop_select_1].PER_WEIGHT)
				#prop_inc_2 = np.sum(sipp_merge[year_select & month_select & prop_select_2].PROP_INC_2 * sipp_merge[month_select & year_select & prop_select_2].PER_WEIGHT)
				# OPTION 2 
				prop_inc_1 = np.sum(sipp_merge[year_select & month_select & prop_select_1].NET_PFT_LOSS_BUS_1 * sipp_merge[month_select & year_select & prop_select_1].PER_WEIGHT)
				prop_inc_2 = np.sum(sipp_merge[year_select & month_select & prop_select_2].NET_PFT_LOSS_BUS_2 * sipp_merge[month_select & year_select & prop_select_2].PER_WEIGHT)


				prop_inc[iy-min_year,im-1] = prop_inc_1 + prop_inc_2 # sum over each year with complete info

		prop_inc_annual 		= np.sum(prop_inc, axis=1)/10**9/4 # sum income across 12 months, divide by 4 bec value is for ref period
		prop_inc_annual[0] 		= 0
		prop_inc_annual[-1] 	= 0

		num_uninc_owner     = partner_num_annual+sole_proprietor_num_annual/10**3
		uninc_inc_per_owner = prop_inc_annual/num_uninc_owner

	num_uninc_owner     = partner_num_annual+sole_proprietor_num_annual/10**3


	print('--- Business income and number of owners (in millions) for unincorporated businesses for this SIPP panel across years ---')
	print("Years %s ---" % np.arange(min_year,max_year+1))
	print("Number of owner is %s ---" % num_uninc_owner)
	if sipp_panel >= 2004:
		print("Business income is %s ---" % prop_inc_annual)
		print("Business income per owner is %s ---" % uninc_inc_per_owner)
	else:
		print("Business income not reported for this panel")
		print("Business income per owner not reported for this panel")



	############ DIVIDEND YIELDS ############
	### Construct annualized business income ###

	if sipp_panel == 2004:
		div_yield_year_mat = [2004,2005]
	elif sipp_panel == 2008:
		div_yield_year_mat = [2009,2010,2011]
	else:
		div_yield_year_mat = []


	for div_yield_year in div_yield_year_mat:
		# Convert to monthly
		sipp_merge['NET_PFT_LOSS_BUS_1_MONTHLY'] = sipp_merge['NET_PFT_LOSS_BUS_1']/4
		sipp_merge['NET_PFT_LOSS_BUS_2_MONTHLY'] = sipp_merge['NET_PFT_LOSS_BUS_2']/4

		sipp_merge['NON_ZERO_BUS_1'] = 0
		sipp_merge.NON_ZERO_BUS_1.loc[ (sipp_merge.NET_PFT_LOSS_BUS_1_MONTHLY!=0)] = 1 
		sipp_merge['NON_ZERO_BUS_2'] = 0
		sipp_merge.NON_ZERO_BUS_2.loc[ (sipp_merge.NET_PFT_LOSS_BUS_2_MONTHLY!=0)] = 1 

		sum_inc = lambda x: x.sum()
		sipp_merge['NET_PFT_LOSS_BUS_1_ANN0'] = (sipp_merge.groupby(['PID','year'])['NET_PFT_LOSS_BUS_1_MONTHLY'].transform(sum_inc)) # do it for the year
		sipp_merge['NET_PFT_LOSS_BUS_2_ANN0'] = (sipp_merge.groupby(['PID','year'])['NET_PFT_LOSS_BUS_2_MONTHLY'].transform(sum_inc)) # do it for the year

		sipp_merge['NON_ZERO_BUS_1_ANN'] = (sipp_merge.groupby(['PID','year'])['NON_ZERO_BUS_1'].transform(sum_inc)) # do it for the year
		sipp_merge['NON_ZERO_BUS_2_ANN'] = (sipp_merge.groupby(['PID','year'])['NON_ZERO_BUS_2'].transform(sum_inc)) # do it for the year

		# Sum and annualize
		sipp_merge['NET_PFT_LOSS_BUS_1_ANN'] = sipp_merge['NET_PFT_LOSS_BUS_1_ANN0'] * (1/(sipp_merge['NON_ZERO_BUS_1_ANN']/12))
		sipp_merge['NET_PFT_LOSS_BUS_2_ANN'] = sipp_merge['NET_PFT_LOSS_BUS_2_ANN0'] * (1/(sipp_merge['NON_ZERO_BUS_2_ANN']/12))

		sipp_merge['NET_PFT_LOSS_BUS_1_ANN'] = sipp_merge['NET_PFT_LOSS_BUS_1_ANN'].fillna(0)
		sipp_merge['NET_PFT_LOSS_BUS_2_ANN'] = sipp_merge['NET_PFT_LOSS_BUS_2_ANN'].fillna(0)

		# Return

		sipp_merge['BUSS_1_RETURN'] = sipp_merge['NET_PFT_LOSS_BUS_1_ANN'] / (sipp_merge['NW_BUS_1']*sipp_merge['PCT_OWN_BUS_1']/100)
		sipp_merge['BUSS_2_RETURN'] = sipp_merge['NET_PFT_LOSS_BUS_2_ANN'] / (sipp_merge['NW_BUS_2']*sipp_merge['PCT_OWN_BUS_2']/100)

		sipp_soleprop_return_buss1 = sipp_merge['BUSS_1_RETURN'].loc[(sipp_merge['LEG_SOLE_PART_1'] == 1) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]
		sipp_soleprop_return_buss2 = sipp_merge['BUSS_2_RETURN'].loc[(sipp_merge['LEG_SOLE_PART_2'] == 1) & (sipp_merge['NW_BUS_2']>0) & (sipp_merge['year'] == div_yield_year)]

		sipp_partnership_return_buss1 = sipp_merge['BUSS_1_RETURN'].loc[(sipp_merge['LEG_SOLE_PART_1'] == 2) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]
		sipp_partnership_return_buss2 = sipp_merge['BUSS_2_RETURN'].loc[(sipp_merge['LEG_SOLE_PART_2'] == 2) & (sipp_merge['NW_BUS_2']>0) & (sipp_merge['year'] == div_yield_year)]

		sipp_uninc_return_buss1 = sipp_merge['BUSS_1_RETURN'].loc[( (sipp_merge['LEG_SOLE_PART_1'] == 1)|(sipp_merge['LEG_SOLE_PART_1'] == 2) ) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]
		sipp_uninc_return_buss2 = sipp_merge['BUSS_1_RETURN'].loc[( (sipp_merge['LEG_SOLE_PART_1'] == 1)|(sipp_merge['LEG_SOLE_PART_1'] == 2) ) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]

		sipp_soleprop_return 		= np.concatenate([sipp_soleprop_return_buss1,sipp_soleprop_return_buss2])
		sipp_partnership_return 	= np.concatenate([sipp_partnership_return_buss1,sipp_partnership_return_buss2])
		sipp_uninc_return 			= np.concatenate([sipp_uninc_return_buss1,sipp_uninc_return_buss2])

		# Net Worth
		sipp_soleprop_nw_buss1 = sipp_merge['NW_BUS_1'].loc[(sipp_merge['LEG_SOLE_PART_1'] == 1) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]
		sipp_soleprop_nw_buss2 = sipp_merge['NW_BUS_2'].loc[(sipp_merge['LEG_SOLE_PART_2'] == 1) & (sipp_merge['NW_BUS_2']>0) & (sipp_merge['year'] == div_yield_year)]

		sipp_partnership_nw_buss1 = sipp_merge['NW_BUS_1'].loc[(sipp_merge['LEG_SOLE_PART_1'] == 2) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]
		sipp_partnership_nw_buss2 = sipp_merge['NW_BUS_2'].loc[(sipp_merge['LEG_SOLE_PART_2'] == 2) & (sipp_merge['NW_BUS_2']>0) & (sipp_merge['year'] == div_yield_year)]

		sipp_uninc_nw_buss1 = sipp_merge['NW_BUS_1'].loc[( (sipp_merge['LEG_SOLE_PART_1'] == 1)|(sipp_merge['LEG_SOLE_PART_1'] == 2) ) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]
		sipp_uninc_nw_buss2 = sipp_merge['NW_BUS_2'].loc[( (sipp_merge['LEG_SOLE_PART_1'] == 1)|(sipp_merge['LEG_SOLE_PART_1'] == 2) ) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]

		sipp_soleprop_nw 		= np.concatenate([sipp_soleprop_nw_buss1,sipp_soleprop_nw_buss2])
		sipp_partnership_nw 	= np.concatenate([sipp_partnership_nw_buss1,sipp_partnership_nw_buss2])
		sipp_uninc_nw 			= np.concatenate([sipp_uninc_nw_buss1,sipp_uninc_nw_buss2])


		# Income
		sipp_soleprop_income_buss1 = sipp_merge['NET_PFT_LOSS_BUS_1_ANN'].loc[(sipp_merge['LEG_SOLE_PART_1'] == 1) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]
		sipp_soleprop_income_buss2 = sipp_merge['NET_PFT_LOSS_BUS_2_ANN'].loc[(sipp_merge['LEG_SOLE_PART_2'] == 1) & (sipp_merge['NW_BUS_2']>0) & (sipp_merge['year'] == div_yield_year)]

		sipp_partnership_income_buss1 = sipp_merge['NET_PFT_LOSS_BUS_1_ANN'].loc[(sipp_merge['LEG_SOLE_PART_1'] == 2) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]
		sipp_partnership_income_buss2 = sipp_merge['NET_PFT_LOSS_BUS_2_ANN'].loc[(sipp_merge['LEG_SOLE_PART_2'] == 2) & (sipp_merge['NW_BUS_2']>0) & (sipp_merge['year'] == div_yield_year)]

		sipp_uninc_income_buss1 = sipp_merge['NET_PFT_LOSS_BUS_1_ANN'].loc[( (sipp_merge['LEG_SOLE_PART_1'] == 1)|(sipp_merge['LEG_SOLE_PART_1'] == 2) ) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]
		sipp_uninc_income_buss2 = sipp_merge['NET_PFT_LOSS_BUS_2_ANN'].loc[( (sipp_merge['LEG_SOLE_PART_1'] == 1)|(sipp_merge['LEG_SOLE_PART_1'] == 2) ) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]

		sipp_soleprop_income 		= np.concatenate([sipp_soleprop_income_buss1,sipp_soleprop_income_buss2])
		sipp_partnership_income 	= np.concatenate([sipp_partnership_income_buss1,sipp_partnership_income_buss2])
		sipp_uninc_income 			= np.concatenate([sipp_uninc_income_buss1,sipp_uninc_income_buss2])

		# Weighted share
		sipp_merge['BUSS_1_WEIGHTED_SHARE'] = sipp_merge['PER_WEIGHT']*sipp_merge['PCT_OWN_BUS_1']/100
		sipp_merge['BUSS_2_WEIGHTED_SHARE'] = sipp_merge['PER_WEIGHT']*sipp_merge['PCT_OWN_BUS_2']/100

		sipp_soleprop_weighted_share_buss1 = sipp_merge['BUSS_1_WEIGHTED_SHARE'].loc[(sipp_merge['LEG_SOLE_PART_1'] == 1) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]
		sipp_soleprop_weighted_share_buss2 = sipp_merge['BUSS_2_WEIGHTED_SHARE'].loc[(sipp_merge['LEG_SOLE_PART_2'] == 1) & (sipp_merge['NW_BUS_2']>0) & (sipp_merge['year'] == div_yield_year)]

		sipp_partnership_weighted_share_buss1 = sipp_merge['BUSS_1_WEIGHTED_SHARE'].loc[(sipp_merge['LEG_SOLE_PART_1'] == 2) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]
		sipp_partnership_weighted_share_buss2 = sipp_merge['BUSS_2_WEIGHTED_SHARE'].loc[(sipp_merge['LEG_SOLE_PART_2'] == 2) & (sipp_merge['NW_BUS_2']>0) & (sipp_merge['year'] == div_yield_year)]

		sipp_uninc_weighted_share_buss1 = sipp_merge['BUSS_1_WEIGHTED_SHARE'].loc[( (sipp_merge['LEG_SOLE_PART_1'] == 1)|(sipp_merge['LEG_SOLE_PART_1'] == 2) ) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]
		sipp_uninc_weighted_share_buss2 = sipp_merge['BUSS_2_WEIGHTED_SHARE'].loc[( (sipp_merge['LEG_SOLE_PART_1'] == 1)|(sipp_merge['LEG_SOLE_PART_1'] == 2) ) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]

		sipp_soleprop_weighted_share 		= np.concatenate([sipp_soleprop_weighted_share_buss1,sipp_soleprop_weighted_share_buss2])
		sipp_partnership_weighted_share 	= np.concatenate([sipp_partnership_weighted_share_buss1,sipp_partnership_weighted_share_buss2])
		sipp_uninc_weighted_share 			= np.concatenate([sipp_uninc_weighted_share_buss1,sipp_uninc_weighted_share_buss2])


		# Weight
		sipp_soleprop_weight_buss1 = sipp_merge['PER_WEIGHT'].loc[(sipp_merge['LEG_SOLE_PART_1'] == 1) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]
		sipp_soleprop_weight_buss2 = sipp_merge['PER_WEIGHT'].loc[(sipp_merge['LEG_SOLE_PART_2'] == 1) & (sipp_merge['NW_BUS_2']>0) & (sipp_merge['year'] == div_yield_year)]

		sipp_partnership_weight_buss1 = sipp_merge['PER_WEIGHT'].loc[(sipp_merge['LEG_SOLE_PART_1'] == 2) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]
		sipp_partnership_weight_buss2 = sipp_merge['PER_WEIGHT'].loc[(sipp_merge['LEG_SOLE_PART_2'] == 2) & (sipp_merge['NW_BUS_2']>0) & (sipp_merge['year'] == div_yield_year)]

		sipp_uninc_weight_buss1 = sipp_merge['PER_WEIGHT'].loc[( (sipp_merge['LEG_SOLE_PART_1'] == 1)|(sipp_merge['LEG_SOLE_PART_1'] == 2) ) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]
		sipp_uninc_weight_buss2 = sipp_merge['PER_WEIGHT'].loc[( (sipp_merge['LEG_SOLE_PART_1'] == 1)|(sipp_merge['LEG_SOLE_PART_1'] == 2) ) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]

		sipp_soleprop_weight 		= np.concatenate([sipp_soleprop_weight_buss1,sipp_soleprop_weight_buss2])
		sipp_partnership_weight 	= np.concatenate([sipp_partnership_weight_buss1,sipp_partnership_weight_buss2])
		sipp_uninc_weight 			= np.concatenate([sipp_uninc_weight_buss1,sipp_uninc_weight_buss2])


		# Share
		sipp_merge['SHARE_OWN_BUS_1'] = sipp_merge['PCT_OWN_BUS_1']/100
		sipp_merge['SHARE_OWN_BUS_2'] = sipp_merge['PCT_OWN_BUS_2']/100

		sipp_soleprop_share_buss1 = sipp_merge['SHARE_OWN_BUS_1'].loc[(sipp_merge['LEG_SOLE_PART_1'] == 1) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]
		sipp_soleprop_share_buss2 = sipp_merge['SHARE_OWN_BUS_2'].loc[(sipp_merge['LEG_SOLE_PART_2'] == 1) & (sipp_merge['NW_BUS_2']>0) & (sipp_merge['year'] == div_yield_year)]

		sipp_partnership_share_buss1 = sipp_merge['SHARE_OWN_BUS_1'].loc[(sipp_merge['LEG_SOLE_PART_1'] == 2) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]
		sipp_partnership_share_buss2 = sipp_merge['SHARE_OWN_BUS_2'].loc[(sipp_merge['LEG_SOLE_PART_2'] == 2) & (sipp_merge['NW_BUS_2']>0) & (sipp_merge['year'] == div_yield_year)]

		sipp_uninc_share_buss1 = sipp_merge['SHARE_OWN_BUS_1'].loc[( (sipp_merge['LEG_SOLE_PART_1'] == 1)|(sipp_merge['LEG_SOLE_PART_1'] == 2) ) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]
		sipp_uninc_share_buss2 = sipp_merge['SHARE_OWN_BUS_2'].loc[( (sipp_merge['LEG_SOLE_PART_1'] == 1)|(sipp_merge['LEG_SOLE_PART_1'] == 2) ) & (sipp_merge['NW_BUS_1']>0) & (sipp_merge['year'] == div_yield_year)]

		sipp_soleprop_share 	= np.concatenate([sipp_soleprop_share_buss1,sipp_soleprop_share_buss2])
		sipp_partnership_share 	= np.concatenate([sipp_partnership_share_buss1,sipp_partnership_share_buss2])
		sipp_uninc_share 		= np.concatenate([sipp_uninc_share_buss1,sipp_uninc_share_buss2])

		# Sole Prop Div Yields
		sipp_soleprop_return_data = pd.DataFrame()
		sipp_soleprop_return_data['sipp_soleprop_return']  			  = sipp_soleprop_return
		sipp_soleprop_return_data['sipp_soleprop_weighted_share'] 	  = sipp_soleprop_weighted_share
		sipp_soleprop_return_data['sipp_soleprop_weight'] 	  		  = sipp_soleprop_weight
		sipp_soleprop_return_data['sipp_soleprop_nw'] 	  			  = sipp_soleprop_nw
		sipp_soleprop_return_data['sipp_soleprop_share'] 	  		  = sipp_soleprop_share
		sipp_soleprop_return_data['sipp_soleprop_buss_nw'] 	  		  = sipp_soleprop_nw * (1/sipp_soleprop_share) # CHECK FOR 0 SHARES BUT STILL OWNERS
		sipp_soleprop_return_data['sipp_soleprop_income'] 	  		  = sipp_soleprop_income 

		temp_soleprop           = wquantile(sipp_soleprop_return_data[['sipp_soleprop_buss_nw','sipp_soleprop_weighted_share']],0.01)
		bottom1_soleprop        = temp_soleprop.iloc[0]
		sipp_soleprop_return_data['high_buss_nw_ind'] = 0
		sipp_soleprop_return_data.high_buss_nw_ind.loc[ (sipp_soleprop_return_data.sipp_soleprop_buss_nw>bottom1_soleprop) ] = 1

		# value-weighted average
		sipp_soleprop_value_weighted_av = 100* ( ( sipp_soleprop_return_data['sipp_soleprop_income']*sipp_soleprop_return_data['sipp_soleprop_weight']*sipp_soleprop_return_data['high_buss_nw_ind'] ).sum() ) / ( ( sipp_soleprop_return_data['sipp_soleprop_nw']*sipp_soleprop_return_data['sipp_soleprop_share']*sipp_soleprop_return_data['sipp_soleprop_weight'] * sipp_soleprop_return_data['high_buss_nw_ind'] ).sum() )

		# equal-weighted average
		sipp_soleprop_equal_weighted_av = 100* ( ( ( sipp_soleprop_return_data['sipp_soleprop_income']*sipp_soleprop_return_data['high_buss_nw_ind'] / (sipp_soleprop_return_data['sipp_soleprop_nw']*sipp_soleprop_return_data['sipp_soleprop_share'])  ) * sipp_soleprop_return_data['sipp_soleprop_weighted_share'] ).sum()  ) / ( (sipp_soleprop_return_data['sipp_soleprop_weighted_share']*sipp_soleprop_return_data['high_buss_nw_ind']).sum()  )

		# distribution
		p10_div_yield_soleprop = wquantile(sipp_soleprop_return_data[['sipp_soleprop_return','sipp_soleprop_weighted_share']].loc[ (sipp_soleprop_return_data.high_buss_nw_ind==1 )] , 0.10) * 100
		p25_div_yield_soleprop = wquantile(sipp_soleprop_return_data[['sipp_soleprop_return','sipp_soleprop_weighted_share']].loc[ (sipp_soleprop_return_data.high_buss_nw_ind==1 )] , 0.25) * 100
		p50_div_yield_soleprop = wquantile(sipp_soleprop_return_data[['sipp_soleprop_return','sipp_soleprop_weighted_share']].loc[ (sipp_soleprop_return_data.high_buss_nw_ind==1 )] , 0.50) * 100
		p75_div_yield_soleprop = wquantile(sipp_soleprop_return_data[['sipp_soleprop_return','sipp_soleprop_weighted_share']].loc[ (sipp_soleprop_return_data.high_buss_nw_ind==1 )] , 0.75) * 100
		p90_div_yield_soleprop = wquantile(sipp_soleprop_return_data[['sipp_soleprop_return','sipp_soleprop_weighted_share']].loc[ (sipp_soleprop_return_data.high_buss_nw_ind==1 )] , 0.90) * 100
		'''
		print('Sole Proprietorship')
		print(sipp_soleprop_value_weighted_av)
		print(sipp_soleprop_equal_weighted_av)
		#print(p10_div_yield_soleprop.iloc[0])
		print(p25_div_yield_soleprop.iloc[0])
		print(p50_div_yield_soleprop.iloc[0])
		print(p75_div_yield_soleprop.iloc[0])
		print(p90_div_yield_soleprop.iloc[0])
		'''
		# Partnerships
		sipp_partnership_return_data = pd.DataFrame()
		sipp_partnership_return_data['sipp_partnership_return']  			  = sipp_partnership_return
		sipp_partnership_return_data['sipp_partnership_weighted_share'] 	  = sipp_partnership_weighted_share
		sipp_partnership_return_data['sipp_partnership_weight'] 	  		  = sipp_partnership_weight
		sipp_partnership_return_data['sipp_partnership_nw'] 	  			  = sipp_partnership_nw
		sipp_partnership_return_data['sipp_partnership_share'] 	  		  	  = sipp_partnership_share
		sipp_partnership_return_data['sipp_partnership_buss_nw'] 	  		  = sipp_partnership_nw * (1/sipp_partnership_share) # CHECK FOR 0 SHARES BUT STILL OWNERS
		sipp_partnership_return_data['sipp_partnership_income'] 	  		  = sipp_partnership_income 

		temp_partnership           = wquantile(sipp_partnership_return_data[['sipp_partnership_buss_nw','sipp_partnership_weighted_share']],0.01)
		bottom1_partnership        = temp_partnership.iloc[0]
		sipp_partnership_return_data['high_buss_nw_ind'] = 0
		sipp_partnership_return_data.high_buss_nw_ind.loc[ (sipp_partnership_return_data.sipp_partnership_buss_nw>bottom1_partnership) ] = 1

		# value-weighted average
		sipp_partnership_value_weighted_av = 100* ( ( sipp_partnership_return_data['sipp_partnership_income']*sipp_partnership_return_data['sipp_partnership_weight']*sipp_partnership_return_data['high_buss_nw_ind'] ).sum() ) / ( ( sipp_partnership_return_data['sipp_partnership_nw']*sipp_partnership_return_data['sipp_partnership_share']*sipp_partnership_return_data['sipp_partnership_weight'] * sipp_partnership_return_data['high_buss_nw_ind'] ).sum() )

		# equal-weighted average
		sipp_partnership_equal_weighted_av = 100* ( ( ( sipp_partnership_return_data['sipp_partnership_income']*sipp_partnership_return_data['high_buss_nw_ind'] / (sipp_partnership_return_data['sipp_partnership_nw']*sipp_partnership_return_data['sipp_partnership_share'])  ) * sipp_partnership_return_data['sipp_partnership_weighted_share'] ).sum()  ) / ( (sipp_partnership_return_data['sipp_partnership_weighted_share']*sipp_partnership_return_data['high_buss_nw_ind']).sum()  )

		# distribution
		p10_div_yield_partnership = wquantile(sipp_partnership_return_data[['sipp_partnership_return','sipp_partnership_weighted_share']].loc[ (sipp_partnership_return_data.high_buss_nw_ind==1 )] , 0.10) * 100
		p25_div_yield_partnership = wquantile(sipp_partnership_return_data[['sipp_partnership_return','sipp_partnership_weighted_share']].loc[ (sipp_partnership_return_data.high_buss_nw_ind==1 )] , 0.25) * 100
		p50_div_yield_partnership = wquantile(sipp_partnership_return_data[['sipp_partnership_return','sipp_partnership_weighted_share']].loc[ (sipp_partnership_return_data.high_buss_nw_ind==1 )] , 0.50) * 100
		p75_div_yield_partnership = wquantile(sipp_partnership_return_data[['sipp_partnership_return','sipp_partnership_weighted_share']].loc[ (sipp_partnership_return_data.high_buss_nw_ind==1 )] , 0.75) * 100
		p90_div_yield_partnership = wquantile(sipp_partnership_return_data[['sipp_partnership_return','sipp_partnership_weighted_share']].loc[ (sipp_partnership_return_data.high_buss_nw_ind==1 )] , 0.90) * 100

		'''
		print('Partnership')
		print(sipp_partnership_value_weighted_av)
		print(sipp_partnership_equal_weighted_av)
		#print(p10_div_yield_partnership.iloc[0])
		print(p25_div_yield_partnership.iloc[0])
		print(p50_div_yield_partnership.iloc[0])
		print(p75_div_yield_partnership.iloc[0])
		print(p90_div_yield_partnership.iloc[0])
		'''

		# Uninc
		sipp_uninc_return_data = pd.DataFrame()
		sipp_uninc_return_data['sipp_uninc_return']  			  = sipp_uninc_return
		sipp_uninc_return_data['sipp_uninc_weighted_share'] 	  = sipp_uninc_weighted_share
		sipp_uninc_return_data['sipp_uninc_weight'] 	  		  = sipp_uninc_weight
		sipp_uninc_return_data['sipp_uninc_nw'] 	  			  = sipp_uninc_nw
		sipp_uninc_return_data['sipp_uninc_share'] 	  		  	  = sipp_uninc_share
		sipp_uninc_return_data['sipp_uninc_buss_nw'] 	  		  = sipp_uninc_nw * (1/sipp_uninc_share) # CHECK FOR 0 SHARES BUT STILL OWNERS
		sipp_uninc_return_data['sipp_uninc_income'] 	  		  = sipp_uninc_income 

		temp_uninc           = wquantile(sipp_uninc_return_data[['sipp_uninc_buss_nw','sipp_uninc_weighted_share']],0.01)
		bottom1_uninc        = temp_uninc.iloc[0]
		sipp_uninc_return_data['high_buss_nw_ind'] = 0
		sipp_uninc_return_data.high_buss_nw_ind.loc[ (sipp_uninc_return_data.sipp_uninc_buss_nw>bottom1_uninc) ] = 1

		# value-weighted average
		sipp_uninc_value_weighted_av = 100* ( ( sipp_uninc_return_data['sipp_uninc_income']*sipp_uninc_return_data['sipp_uninc_weight']*sipp_uninc_return_data['high_buss_nw_ind'] ).sum() ) / ( ( sipp_uninc_return_data['sipp_uninc_nw']*sipp_uninc_return_data['sipp_uninc_share']*sipp_uninc_return_data['sipp_uninc_weight'] * sipp_uninc_return_data['high_buss_nw_ind'] ).sum() )

		# equal-weighted average
		sipp_uninc_equal_weighted_av = 100* ( ( ( sipp_uninc_return_data['sipp_uninc_income']*sipp_uninc_return_data['high_buss_nw_ind'] / (sipp_uninc_return_data['sipp_uninc_nw']*sipp_uninc_return_data['sipp_uninc_share'])  ) * sipp_uninc_return_data['sipp_uninc_weighted_share'] ).sum()  ) / ( (sipp_uninc_return_data['sipp_uninc_weighted_share']*sipp_uninc_return_data['high_buss_nw_ind']).sum()  )

		# distribution
		p10_div_yield_uninc = wquantile(sipp_uninc_return_data[['sipp_uninc_return','sipp_uninc_weighted_share']].loc[ (sipp_uninc_return_data.high_buss_nw_ind==1 )] , 0.10) * 100
		p25_div_yield_uninc = wquantile(sipp_uninc_return_data[['sipp_uninc_return','sipp_uninc_weighted_share']].loc[ (sipp_uninc_return_data.high_buss_nw_ind==1 )] , 0.25) * 100
		p50_div_yield_uninc = wquantile(sipp_uninc_return_data[['sipp_uninc_return','sipp_uninc_weighted_share']].loc[ (sipp_uninc_return_data.high_buss_nw_ind==1 )] , 0.50) * 100
		p75_div_yield_uninc = wquantile(sipp_uninc_return_data[['sipp_uninc_return','sipp_uninc_weighted_share']].loc[ (sipp_uninc_return_data.high_buss_nw_ind==1 )] , 0.75) * 100
		p90_div_yield_uninc = wquantile(sipp_uninc_return_data[['sipp_uninc_return','sipp_uninc_weighted_share']].loc[ (sipp_uninc_return_data.high_buss_nw_ind==1 )] , 0.90) * 100


		print("--- Net income yields for unincorporated businesses for year %s ---" % div_yield_year)
		print("Value weighted average is %s ---" % sipp_uninc_value_weighted_av)
		print("Equally weighted average is %s ---" % sipp_uninc_equal_weighted_av)
		print("p25 of the net income yield distribution is %s ---" % p25_div_yield_uninc.iloc[0])
		print("p50 of the net income yield distribution is %s ---" % p50_div_yield_uninc.iloc[0])
		print("p75 of the net income yield distribution is %s ---" % p75_div_yield_uninc.iloc[0])


		my_dataframe = pd.read_csv(sipp_dir+'ta11.csv', index_col=0)
		my_dataframe.loc[div_yield_year] = [sipp_soleprop_value_weighted_av, sipp_soleprop_equal_weighted_av, p25_div_yield_soleprop.iloc[0], p50_div_yield_soleprop.iloc[0], p75_div_yield_soleprop.iloc[0],
		sipp_partnership_value_weighted_av, sipp_partnership_equal_weighted_av, p25_div_yield_partnership.iloc[0], p50_div_yield_partnership.iloc[0], p75_div_yield_partnership.iloc[0],
		sipp_uninc_value_weighted_av, sipp_uninc_equal_weighted_av, p25_div_yield_uninc.iloc[0], p50_div_yield_uninc.iloc[0], p75_div_yield_uninc.iloc[0]]
		
		if div_yield_year == 2011:
			my_dataframe.loc["Mean"] = my_dataframe.mean()

		my_dataframe.to_csv(sipp_dir+'ta11.csv', index = True)