Many (but not all) economists love the market as an organizing institution. This affection for markets is traced back by some to Adam Smith and his "invisible hand". Others may not think that the market is so great, but basically distrust other organizing institutions like a government. They will point, for instance, to the fall of communism in the former Soviet Union.

Although you, as an economist, have seen many market models already, it may be illustrative to go back to basics. To understand the advantages of the market, let us consider a very simple economy. We focus on one type of product and there is an exogenous endowment of this product (supply) equal to number_of_goods. Further, in this economy there are number_of_agents agents who have a valuation for this good which is randomly drawn from a normal distribution. An agent can at max. consume one unit of the product and her utility is then given by her valuation. The value of consuming zero units and the additional value of consuming more than one unit both equal 0.

The vector valuations contains for each agent her valuation and we sort this vector such that the first agent has the highest valuation.

number_of_agents = 1000
number_of_goods = 100

valuations = sorted(pm.Normal.dist(100,20).random(size=number_of_agents),reverse = True)

Note that we are using pymc3 here to generate random numbers from a distribution. There are also other python libraries which can do this, e.g scipy. We use pymc3 as we will use it later for estimation as well.

Question What is the economic name for the following expression? To answer this question, you need to understand both how indexing works in python and which economic concept is captured by this expression.

print("{0:.2f}".format(valuations[number_of_goods]))

125.21

Suppose that we have an omniscient social planner who knows the valuation of every agent in the economy. This planner aims to maximize the (unweighted) sum of agents' utilities.

Question Calculate the total welfare that this planner can achieve. Denote this value max_welfare.

Hence, this is the best that we can do. It gives us an upperbound on the welfare that can be achieved. Having an omniscient social planner seems unrealistic, but perhaps there is an institution that can achieve this outcome without omniscient intervention. You guessed it…

Although the results above look great, the assumptions we made, may not be realistic in every market. To illustrate, without saying so, we assumed above that the market is perfectly competitive and without external effects. Here we program three reasons why the market outcome may not necessarily lead to maximum welfare. First, we look at income inequality and the problem that this causes for the market. Then we consider market power and finally we model external effects.

Question What is adverse selection?

Consider an economy with number_of_agents agents. Each agent has an endowment/income equal to income and faces a potential loss of the size cost. Agents differ in the probability \(\pi\) of this loss. We randomly draw 50 values for \(\pi\) assuming it is uniformly distributed on \([0,1]\).

Further, agents have a utility function of the form \(u(x)=x^{\rho}\).

income = 1.1
cost = 1
ρ = 0.1
def u(x):
    return x**ρ

number_of_agents = 50

π = pm.Uniform.dist(0.0,1.0).random(size = number_of_agents)
π.sort()

Since we assume that \(\rho \in \langle 0, 1 \rangle\), agents are risk averse and would like to buy insurance which covers the loss. We assume that insurance covers the loss completely at a premium \(\sigma\). As we assume that the probability of loss, \(\pi\), is exogenous, there is no reason to have co-payments of any sort.

An agent buys insurance if and only if

Question Define a function insurance_demand that returns the number of agents buying insurance as a function of the premium \(\sigma\).

We assume that this insurance market is perfectly competitive. That is, for each quantity supplied, insurance companies compete down the price such that the premium equals the average (expected) cost of the agents buying insurance.

Question Explain the code of the following function.

def insurance_supply(Q):
    return np.mean(π[-Q:])*cost

We plot demand and supply in one figure. In addition, we plot the marginal costs curve.

range_Q = np.arange(1,number_of_agents+1,1)
range_sigma = np.arange(0,1.01,0.01)
plt.plot(range_Q,[insurance_supply(Q) for Q in range_Q],label="insurance supply")
plt.plot([insurance_demand(sigma) for sigma in range_sigma],range_sigma,label="insurance demand")
plt.plot(range_Q,[π[-Q]*cost for Q in range_Q],label="marginal cost")
plt.legend()
plt.xlabel('$Q$')
plt.ylabel('$\sigma$')
plt.title('Perfectly competitive insurance market')
plt.show()

<Figure size 432x288 with 1 Axes>

Question Interpret this figure. In particular,

• explain why all curves are downward sloping (is supply not usually upward sloping?)
• what is approx. the equilibrium premium \(\sigma\)?
• is the market outcome efficient?
• what can we learn from the marginal cost curve?

Assignment Show graphically the effect of an increase in income on the market outcome. Does the inefficiency increase or decrease with income? Why?

With moral hazard, agents take hidden actions. The actions that they take are affected by the incentives that they face. We consider this in the context of taxation.

People differ in their productivity. For some people it is easy to generate a gross income \(x\), for others generating such an income would be very costly in terms of effort. In the real world, such differences in productivity can be caused by IQ, education, health status etc. Here, we simply model this as an effort cost. People with a high effort cost have lower productivity than people with low effort costs. We assume that the effort cost is log-normally distributed.

The government uses a linear tax schedule: \(\tau x - \tau_0\). Hence, when you have a gross income \(x\), your net income equals \((1-\tau)x+\tau_0\). Where we assume that for the economy as a whole the tax revenue is redistributed among the population. Hence, number_of_agents times \(\tau_0\) has to equal the total revenue from the marginal tax rate \(\tau\).

Agents maximize their utility by choosing production \(x\):

where agents differ in \(c\) and \(c\) is not observable.

These two aspects are important: if \(c\) were observable or if everyone was symmetric (had the same \(c\)) taxation would be easy. To see why, first note that income \(x\) is apparently observable since taxation depends on it. Hence, the government could say to an agent \(c\): I want you to produce income \(x\) and you give me a share \(\tau\) of this income.

In our set-up with heterogeneity in \(c\) and \(c\) unobservable, the government cannot force people to generate income \(x\) because some of these agents may have such a high \(c\) that this is inefficient (or even impossible).

Hence, the government sets the tax schedule (in our case linear) and allows each agent to choose her own production level. The higher \(\tau\), the lower an agent's production will be.

number_of_agents = 200
effort_costs = pm.Lognormal.dist(mu=0.0,sd=0.5).random(size=number_of_agents)
def effort(c,τ):
    sol = optimize.minimize(lambda x: -(x*(1-τ)-c*x**2),1)
    return sol.x

We use the following welfare function:

With \(\rho=1\), the social planner just maximizes the sum of utility. With \(\rho<1\), the planner has a taste for redistribution: agents with low utility get a relatively high weight in this welfare function.

The function Welfare first calculates for a given \(\tau\), what the value of \(\tau_0\) is (using budget balance for the government). Then for these values of \(\tau\) and \(\tau_{0}\), \(W\) is calculated.

def Welfare(τ,ρ):
    τ_0 = np.mean([τ*effort(c,τ) for c in effort_costs])
    return (np.sum([((1-τ)*effort(c,τ)+τ_0 - c*effort(c,τ)**2)**ρ for c in effort_costs]))**(1/ρ)

Question Plot Welfare as a function of \(\tau\) for \(\rho=1\). What is the welfare maximizing tax rate? Why?

Question What happens to the optimal tax rate as \(\rho<1\) falls?

Assignment Redefine the function Welfare above such that it uses Rawls' criterion of maximizing the utility of the person who is worse off in society. Further, suppose that the government needs \(g\) per head to finance a public good. What is the effect of \(g\) on the optimal marginal tax rate?

The first problem in financial (and other) markets is limited liability. When banks go bankrupt they "only" loose their equity even if the debts that they accumulated exceed their equity. We run some simulations to show that this leads to banks taking excessive risks from a social point of view.

Let \(x\) denote an investment opportunity: \(x\) is a vector with dimension 1000. That is, we assume that there are 1000 states of the world and \(x\) gives us the return in each of these states of the world. To find the expected profit, we take the average over the return in all these states of the world. However, if \(x\) is "very negative" (a big loss), the bank goes bankrupt and the owners only loose their equity.

def profit(x,equity=0):
    return np.mean(np.maximum(x,-equity))

Question Plot the function profit for the case where equity =10 and \(x\) is a scalar.

Note that although we think of \(x\) as a vector (and we will use this in a second), python does not know this nor does it care. Hence, we can plot profit as a function of the scalar \(x\).

Question Create an investment opportunity vector_returns where the returns are normally distributed with mean \(-10\) and standard deviation 100. As mentioned above, we assume that there are 1000 states of the world.

Question Calculate the expected (social) return from this investment vector_returns.

Question Calculate the expected profits of the vector_returns. Compare the outcome to the one above. What is the interpretation?

Question Explain what is coded in the following code cell. E.g. what is the point of the 10000?

v_std = np.arange(0,200,1)
v_returns = [pm.Normal.dist(-10,std).random(size=1000) for std in v_std]
plt.scatter([np.std(vx) for vx in v_returns],[profit(vx) for vx in v_returns], label="no equity")
plt.scatter([np.std(vx) for vx in v_returns],[profit(vx,60) for vx in v_returns], label="equity equals 60")
plt.scatter([np.std(vx) for vx in v_returns],[profit(vx,10000) for vx in v_returns], label="social value")
plt.xlabel('$\sigma$')
plt.ylabel('return')
plt.legend()
plt.show()

<Figure size 432x288 with 1 Axes>

Question Explain the economic intuition/interpretation of the graph above.

Bonus contracts are generally a reaction to asymmetric information. Let's first consider moral hazard.

A good reason to use python for data analysis is the option to get on-line data directly into your notebook without going to the website first to download this data. A number of institutes have such python API's. To illustrate this, we use the Worldbank API as described on this website.

The advantage of doing your analysis in this way is that your research becomes better reproducible. Everyone can run the same code and then go through your code of data cleaning steps to end up with the same data set. If instead you first download the data to your computer, then use excel to clean the data and then start analyzing it (say, with stata), no one will be able to reproduce the exact steps that you have taken.

To illustrate the Worldbank API, we will look at the development over time of inequality in gdp per head. So we want measures of gdp per head. The API allows us to search for such indicators in the Worldbank data set. The column on the left gives the name of the variables (that we will use below to download the data); the column on the right explains what the variable provides.

wb.search_indicators("gdp per capita")

6.0.GDPpc_constant      	GDP per capita, PPP (constant 2011 international $) 
FB.DPT.INSU.PC.ZS       	Deposit insurance coverage (% of GDP per capita)
NV.AGR.PCAP.KD.ZG       	Real agricultural GDP per capita growth rate (%)
NY.GDP.PCAP.PP.KD.ZG    	GDP per capita, PPP annual growth (%)
NY.GDP.PCAP.PP.KD.87    	GDP per capita, PPP (constant 1987 international $)
NY.GDP.PCAP.PP.KD       	GDP per capita, PPP (constant 2011 international $)
NY.GDP.PCAP.PP.CD       	GDP per capita, PPP (current international $)
NY.GDP.PCAP.KN          	GDP per capita (constant LCU)
NY.GDP.PCAP.KD.ZG       	GDP per capita growth (annual %)
NY.GDP.PCAP.KD          	GDP per capita (constant 2010 US$)
NY.GDP.PCAP.CN          	GDP per capita (current LCU)
NY.GDP.PCAP.CD          	GDP per capita (current US$)
SE.XPD.TERT.PC.ZS       	Government expenditure per student, tertiary (% of GDP per capita)
SE.XPD.SECO.PC.ZS       	Government expenditure per student, secondary (% of GDP per capita)
SE.XPD.PRIM.PC.ZS       	Government expenditure per student, primary (% of GDP per capita)
UIS.XUNIT.GDPCAP.4.FSGOV	Government expenditure per post-secondary non-tertiary student as % of GDP per capita (%)
UIS.XUNIT.GDPCAP.3.FSGOV	Government expenditure per upper secondary student as % of GDP per capita (%)
UIS.XUNIT.GDPCAP.2.FSGOV	Government expenditure per lower secondary student as % of GDP per capita (%)

Let's say that we are interested in "GDP per capita, PPP (constant 2011 international $)", we specify this indicator in a dictionary where the key is the "official name" of the variable and the value is the way that we want to refer to the variable (in this case: "GDP_per_head").

With get_dataframe we actually download the data into the dataframe df_wb. We reset the index in this case (just see what happens to the dataframe if you don't do this). And we look at the first 5 rows to get an idea of what the data are.

indicators = {"NY.GDP.PCAP.PP.KD": "GDP_per_head"}
df_wb = wb.get_dataframe(indicators, convert_date=True)
df_wb.reset_index(inplace = True)
df_wb.head()

      country       date  GDP_per_head
0  Arab World 2017-01-01  15413.791998
1  Arab World 2016-01-01  15500.530523
2  Arab World 2015-01-01  15342.766482
3  Arab World 2014-01-01  15199.008915
4  Arab World 2013-01-01  15174.101703

Question What do the last 10 rows look like?

If you like the dataframe that you have downloaded from the web, you can save it with pandas to_csv. We save the data in the subdirectory "data".

df_wb.to_csv('data/worldbank_data_gdp_per_capita.csv')

Let's compare the distribution of gdp per head in 1990 with the distribution in 2017. In order to illustrate how we can combine dataframes in pandas, we first define separate dataframes for the years 1990 and 2017.

df_1990=df_wb[df_wb['date']=='1990-01-01']
df_2017=df_wb[df_wb['date']=='2017-01-01']

Question What does the dataframe df_1990 look like?

Both dataframes have a column country. Hence, we can merge the dataframes on this column. There are a number of how methods, here we use 'inner' which means that only countries that are present in both datasets will be in df_merged. To distinguish the columns, like GDP_per_head from the two dataframes, we can provide suffixes. All columns from df_1990 (except for country) will be suffixed with '_1990'; and similarly for 2017.

df_merged = pd.merge(df_1990, df_2017, on=['country'], suffixes=['_1990', '_2017'], how='inner')

Question To see how the suffixes work, check what df_merged looks like.

Question Plot GPD per head in 1990 against GDP per head in 2017. What do you conclude about the development in inequality in these 27 years?

You may wonder which observations ("dots") correspond to which countries. For this we need a plotting library that is more sophisticated on interactions than matplotlib. A number of these libraries are available; here we consider bokeh. If you want to know more about bokeh, there is a datacamp course.

from bokeh.io import output_file, show, output_notebook
from bokeh.plotting import figure
from bokeh.models import HoverTool
output_notebook()

hover = HoverTool(tooltips=[
     ('country','@country'),
     ])

plot = figure(tools=[hover])
plot.circle('GDP_per_head_1990','GDP_per_head_2017',
    size=10, source=df_merged)
output_file('inequality.html')
show(plot)

./inequality.html

In the last chapter of the Datacamp course Intermediate Python for Data Science you have seen hacker statistics. The idea is to simulate a random process, say, 10,000 times. Then if you want to know the probability that a certain condition is satisfied (say, number of heads bigger than 10 when you throw a coin 15 times), you calculate the number of outcomes where this condition is satisfied and divide it by 10,000.

If you can program, you can recap all the statistics that you were taught (and probably forgot).

Consider the following simple theoretical framework. People get offered at max. one treatment per year. They decide whether either to accept this treatment or to go without treatment.

The figure below plots costs \(c\) of treatment vs. values \(v\) of treatments. Each treatment is a point in this figure; a combination of \(c\) and \(v\).

The red/blue line is the out-of-pocket payment by an agent facing two deductible levels: 365 and 170 euro resp. Treatments above the red line are always accepted. The value exceeds the out-of-pocket payment for both deductibles. Treatments below the blue line are always rejected: even with the low deductible, the value is below the out-of-pocket payment. Treatments in the yellow area are accepted with the low deductible but are rejected with the high deductible. Hence, to quantify the effect of an increase in the deductible, we want to know the probability that treatments fall in the yellow area. The more treatments in the yellow area, the bigger the fall in health care expenditure in response to an increase in \(d\).

def deductible(c,d):
   return min(c,d)

range_c = np.arange(0,500,0.1)
range_v170 = [deductible(c,170) for c in range_c]
range_v365 = [deductible(c,365) for c in range_c]

plt.plot(range_c,range_v365,'-', color = 'r', linewidth = 2, label = '$d=365$')
plt.plot(range_c,range_v170,'-', color = 'b', linewidth = 2, label = '$d=170$')
plt.legend()
plt.fill_between(range_c, range_v170, range_v365, facecolor='yellow')
plt.ylim(0,500)
plt.xlabel('Cost')
plt.ylabel('Value')
plt.show()

<Figure size 432x288 with 1 Axes>

To find the effect of an increase in deductible, we compare health care expenditures in the years 2011 (deductible was 170 euro) and 2014 (deductible was 365 euro). We use data from Vektis. Download from this website the 'csv' files for 2011 and 2014. To use the code below, download these csv-files into the sub-directory "data" (i.e. "data" is sub-directory of the directory in which this notebook resides).

When you open the csv files, you can see that it uses ";" as separator between columns. Hence, we use pandas' read_csv statement where we specify the separator as ';'. The data contain total cost per postal code area for a number of cost categories. The expenditures are grouped by sex and age.

The function get_data_into_shape does a number of things:

• not all health care cost categories in the data "count" as far as the deductible is concerned. Hence, we select the ones that fall under the deductible and sum these as the relevant total expenditure under the deductible.
• the cost categories are in Dutch, hence we translate the labels into English
• we drop variables that we do not need for the analysis here
• we calculate cost per head per postal code area
• we also introduce the log of health care costs per head
• we turn the variable sex into a category with two values ('M' for males, 'V' for females)
• we drop the age category '91+' and turn the remaining ages into integers
• finally, the function returns this new dataframe.

The function illustrates the data manipulation you can do with pandas. We use the same function for the 2011 and 2014 data, which means that we ignore some columns in the 2014 data that are not available in the 2011 data.

df_2014 = pd.read_csv('data/Vektis Open Databestand Zorgverzekeringswet 2014 - postcode3.csv', sep = ';')

cost_categories_under_deductible = ['KOSTEN_MEDISCH_SPECIALISTISCHE_ZORG', 'KOSTEN_MONDZORG', 'KOSTEN_FARMACIE', 'KOSTEN_HULPMIDDELEN', 'KOSTEN_PARAMEDISCHE_ZORG_FYSIOTHERAPIE', 'KOSTEN_PARAMEDISCHE_ZORG_OVERIG', 'KOSTEN_ZIEKENVERVOER_ZITTEND', 'KOSTEN_ZIEKENVERVOER_LIGGEND', 'KOSTEN_GRENSOVERSCHRIJDENDE_ZORG', 'KOSTEN_OVERIG']

def get_data_into_shape(df):
    df['health_expenditure_under_deductible'] = df[cost_categories_under_deductible].sum(axis=1)
    df = df.rename({
        'GESLACHT':'sex',
        'LEEFTIJDSKLASSE':'age',
        'GEMEENTENAAM':'MUNICIPALITY',
        'AANTAL_BSN':'number_citizens',
        'KOSTEN_MEDISCH_SPECIALISTISCHE_ZORG':'hospital_care',
        'KOSTEN_FARMACIE':'pharmaceuticals',
        'KOSTEN_TWEEDELIJNS_GGZ':'mental_care',
        'KOSTEN_HUISARTS_INSCHRIJFTARIEF':'GP_capitation',
        'KOSTEN_HUISARTS_CONSULT':'GP_fee_for_service',
        'KOSTEN_HUISARTS_OVERIG':'GP_other',
        'KOSTEN_MONDZORG':'dental care',
        'KOSTEN_PARAMEDISCHE_ZORG_FYSIOTHERAPIE':'physiotherapy',
        'KOSTEN_KRAAMZORG':'maternity_care',
        'KOSTEN_VERLOSKUNDIGE_ZORG':'obstetrics'
    }, axis='columns')
    df.drop(['AANTAL_VERZEKERDEJAREN',
             'KOSTEN_HULPMIDDELEN',
             'KOSTEN_PARAMEDISCHE_ZORG_OVERIG',
             'KOSTEN_ZIEKENVERVOER_ZITTEND',
             'KOSTEN_ZIEKENVERVOER_LIGGEND',
             'KOSTEN_GRENSOVERSCHRIJDENDE_ZORG',
             'KOSTEN_OVERIG',
             'KOSTEN_EERSTELIJNS_ONDERSTEUNING'],inplace=True,axis=1)
    df.drop(df.index[[0]], inplace=True)
    df['sex'] = df['sex'].astype('category')
    df['age'] = df['age'].astype('category')
    df['costs_per_head']=df['health_expenditure_under_deductible']/df['number_citizens']
    df['log_costs_per_head']=np.log(1+df['health_expenditure_under_deductible']/df['number_citizens'])
    df = df[(df['age'] != '90+')]
    df['age'] = df['age'].astype(int)
    return df

df_2014 = get_data_into_shape(df_2014)
df_2014.head()

We create costs_per_sex_age which contains the average health care expenditure (averaged over postal code areas) for each combination of sex and age in the data. For this we use pandas method groupby.

The outcome of this groupby we plot below.

costs_per_sex_age = df_2014.groupby(['sex','age'])['costs_per_head'].mean()

We plot the distribution of health care expenditure per head with age for males and females.

fig = plt.figure()
ax = costs_per_sex_age['M'].plot()
ax = costs_per_sex_age['V'].plot()
ax.set_xlabel('age')
ax.set_ylabel('costs per head')
ax.set_title('average costs per age and sex')
ax.legend(['male','female'])

<Figure size 432x288 with 1 Axes>

Question Can you interpret how these costs evolve with age and sex?

Question How can a graph like this help us to determine the effect of \(d\) on health care expenditure?

Above we used pymc3 to generate vectors of productivities, valuations, incomes etc. using probability distributions. Here we go the "other way around". We have here distributions of health care expenditures per head and we want to identify the distributions where these come from. To illustrate this, consider the distribution of (average) costs for 30 year old males. Since, health care costs have a skewed distribution, we actually plot the distribution of log costs.

Question Plot health care cost distributions for different age and sex categories.

df_2014.query('sex=="M" & age=="30"')['log_costs_per_head'].hist(bins=50)

<Figure size 432x288 with 1 Axes>

This distribution looks (sort of) normal. Hence, we assume that for each age and sex category log_costs_per_head are normally distributed. This implies that costs_per_head have a log-normal distribution.

Question Why are there no zeroes in this figure? Lots of people have no health care expenditure in a year? How does our log_costs_per_head deal with zero costs? Should there be a spike at 0?

We focus here on health care costs for women. Clearly, a similar analysis can be done for men. In fact, it is also possible to combine men and women into one analysis with gender fixed effects.

Here we focus on women and introduce age-fixed effects. We assume that observed costs \(z\) are log_costs_per_head which are normally distributed with a mean \(\mu\) and standard deviation \(\sigma\) which both vary with age. We do not know these means and standard deviations μ[age], σ[age] but assume they are drawn from prior distributions. A normal distribution for \(\mu\) and a half-normal distribution for \(\sigma\).

This is called a Bayesian analysis which you probably never saw before. Do not worry about this; you do not need to understand the details of this analysis. If you find it interesting, there are some references at the end.

log_costs_per_age_female = df_2014[df_2014['sex']=='V'].groupby(['age'])['log_costs_per_head'].mean()

log_costs_per_head = df_2014[df_2014['sex']=='V'].log_costs_per_head.values
age = df_2014[df_2014['sex']=='V'].age.values


with pm.Model() as model:
    
    μ = pm.Normal('μ', 8, 3, shape=len(set(age)))
    σ = pm.HalfNormal('σ', 4, shape=len(set(age)))
    z = pm.Normal('z', μ[age], σ[age], observed=log_costs_per_head)

with model:
    trace = pm.sample(4000,step = pm.Metropolis(),start = pm.find_MAP())

summary = pm.summary(trace, varnames=['μ'])

pm.plot_posterior(trace, varnames=['μ'],ref_val = log_costs_per_age_female.values)[0]

The figures above compare the estimated distribution for \(\mu\) for each age category with the observed average expenditure for this age category in the data. The observed average expenditure falls in the middle of this distribution for each age category. This suggests that the estimation makes some sense (although more tests should be done).

The "fun" of Bayesian analysis is that we have a (posterior) distribution of the population parameters \(\mu\). With "classic" econometrics this is not possible, since \(\mu\) is not random in that analysis.

To elaborate a little, in classical econometric analysis we consider the probability \(p(\text{data}|\theta)\) where \(\theta\) denotes the parameters of the underlying distribution. Hence, we do \(p(m>10|\mu=9)\). With Bayes, we do \(p(\theta|\text{data})\). The idea is that the data is fixed (we have observed it) but we face uncertainty over the values of the parameters \(\theta\). The graphs above plot for each age the distribution of \(\mu\).

Now we can plot the average \(\mu\) for each age and the observed average expenditure per age category in a graph:

plt.plot(summary['mean'].values,label='calculated means')
plt.plot(log_costs_per_age_female,'o',label='observed means')
plt.legend()

<Figure size 432x288 with 1 Axes>

To see the effect of the deductible, we compare the average \(\mu\) for 17 year olds with the average \(\mu\) for 19 year olds:

summary['mean']['μ__17'] - summary['mean']['μ__19']

0.29334615011390586

This is positive: 17 year olds spend more on health care than 19 year olds. Now you may think that although the average \(\mu\) is higher for 17 year olds than for 19 year olds, but perhaps it is still likely that \(\mu_{19}>\mu_{17}\). Hence we plot these \(\mu\) distributions:

plt.hist(trace['μ'][:,17],density=True,label='age 17',bins=50)
plt.hist(trace['μ'][:,19],density=True,label='age 19',bins=50)
plt.legend()

<Figure size 432x288 with 1 Axes>

Question Does this mean that expenditure of every 17 year old exceeds the expenditure of every 19 year old?

Questions Focus on 19 year olds. What is the probability that \(\mu_{19}>6.45\)?

One explanation for the figure is that 17 year olds do not face a deductible, while 19 year olds do. But an other explanation is that health care expenditure simply differs by age (irrespective of a deductible). In order to control for the age effect, we re-do the analysis above for 2011. Also in 2011 we can take the difference in means for 17 and 19 year olds. If there is (only) a biological reason for the different expenditures between 17 and 19 year olds, the difference in 2011 should be the same as the difference in 2014.

If, however, the difference in expenditures is caused by the deductible, we expect a bigger difference in 2014 than in 2011 as the deductible was higher in 2014 than in 2011. In terms of our model in section 7.1 above: the yellow area is the additional effect due to the higher deductible in 2014.

Hence, we do the same analysis as above for 2011.

df_2011 = pd.read_csv('data/Vektis Open Databestand Zorgverzekeringswet 2011 - postcode3.csv', sep = ';')

df_2011 = get_data_into_shape(df_2011)
df_2011.head()

As above we estimate the model.

log_costs_per_age_female = df_2011[df_2011['sex']=='V'].groupby(['age'])['log_costs_per_head'].mean()

log_costs_per_head = df_2011[df_2011['sex']=='V'].log_costs_per_head.values
age = df_2011[df_2011['sex']=='V'].age.values


with pm.Model() as model_2011:
    
    μ = pm.Normal('μ', 8, 3, shape=len(set(age)))
    σ = pm.HalfCauchy('σ', 4, shape=len(set(age)))
    z = pm.Normal('z', μ[age], σ[age], observed=log_costs_per_head)

with model_2011:
    trace_2011 = pm.sample(4000,step = pm.Metropolis(),start = pm.find_MAP())

summary_2011 = pm.summary(trace_2011, varnames=['μ'])

Recall that in 2014 the difference in mean log expenditure equals: 0.29

summary['mean']['μ__17']-summary['mean']['μ__19']

0.29334615011390586

In 2011, this difference is smaller: 0.19

summary_2011['mean']['μ__17']-summary_2011['mean']['μ__19']

0.19337257026002064

Also this we can do with distributions.

plt.hist(trace['μ'][:,17]-trace['μ'][:,19],density=True,label='difference in 2014',bins=50)
plt.hist(trace_2011['μ'][:,17]-trace_2011['μ'][:,19],density=True,label='difference in 2011',bins=50)

plt.legend()

<Figure size 432x288 with 1 Axes>

Question Calculate the probability that in fact the parameter difference \(\mu_{17}-\mu_{19}\) is bigger in 2011 than in 2014.

Hence we find that government regulation in terms of a deductible has an effect on health care expenditures. People spend less on health care if they face a higher deductible.

One can extend the analysis above to make it more convincing. Think of things like

• add more years (using the Vektis website)
• add year dummies to distinguish year effects (like changes in treatments covered by basic insurance) from changes in deductible
• add the level of the deductible in the estimation.

If you like this way of modelling with pymc3, you can look at this video

A simple introduction to Bayesian estimation with python can be found here. Allen Downey also has free books on statistics with python.

Bibliography

• [tirole_2017] Jean Tirole, Economics for the Common Good, Princeton University Press (2017).
• [bijlsma2018] Bijlsma, Boone & Zwart, Competition for traders and risk, RAND Journal of Economics, 34(4), 737-763 (forthcoming).