Full hierarchical model

foo.py

# For use in Jupyter Notebook include the next line
%matplotlib inline
import pymc3 as pm

degree_indexes = degree_index['index'].values
degree_count = len(degree_indexes)
degree_state_indexes = degree_state_indexes_df['index_d'].values
degree_state_count = len(degree_state_indexes)
degree_state_county_indexes = degree_state_county_indexes_df['index_ds'].values
degree_state_county_count = len(degree_state_county_indexes)

with pm.Model() as model:
    global_m = pm.Normal('global_m', mu=0, sd=100**2)
    global_m_sd = pm.Uniform('global_m_sd', lower=0, upper=1000)
    global_b = pm.Normal('global_b', mu=0, sd=100**2)
    global_b_sd = pm.Uniform('global_b_sd', lower=0, upper=1000)
    
    degree_m = pm.Normal('degree_m', mu=global_m, sd=global_m_sd, shape=degree_count)
    degree_m_sd = pm.Uniform('degree_m_sd', lower=0, upper=1000, shape=degree_count)
    degree_b = pm.Normal('degree_b', mu=global_b, sd=global_b_sd, shape=degree_count)
    degree_b_sd = pm.Uniform('degree_b_sd', lower=0, upper=1000, shape=degree_count)
    
    degree_state_m = pm.Normal('degree_state_m', mu=degree_m[degree_state_indexes], sd=degree_m_sd[degree_state_indexes], shape=degree_state_count)
    degree_state_m_sd = pm.Uniform('degree_state_m_sd', lower=0, upper=1000, shape=degree_state_count)
    degree_state_b = pm.Normal('degree_state_b', mu=degree_b[degree_state_indexes], sd=degree_b_sd[degree_state_indexes], shape=degree_state_count)
    degree_state_b_sd = pm.Uniform('degree_state_b_sd', lower=0, upper=1000, shape=degree_state_count)
    
    degree_state_county_m = pm.Normal('degree_state_county_m', mu=degree_state_m[degree_state_county_indexes], sd=degree_state_m_sd[degree_state_county_indexes], shape=degree_state_county_count)
    degree_state_county_b = pm.Normal('degree_state_county_b', mu=degree_state_b[degree_state_county_indexes], sd=degree_state_b_sd[degree_state_county_indexes], shape=degree_state_county_count)
    
    error = pm.Uniform('error', lower=0, upper=10000)
    
    y_prediction = degree_state_county_m[indexed_salary_df['index'].values] * indexed_salary_df.month_index.values + degree_state_county_b[indexed_salary_df['index'].values]
    
    pm.StudentT('y_like', nu=1, mu=0, sd=error, observed=y_prediction - indexed_salary_df.salary.values)
    
    model_trace = pm.sample(5000)
    
pm.traceplot(model_trace[1000:]);