data = ps.pdio.read_files(ps.examples.get_path('NAT.shp'))
W = ps.queen_from_shapefile(ps.examples.get_path('NAT.shp'))
W.transform = 'r'

data.head()

"""
Explanation: A well-used functionality in PySAL is the use of PySAL to conduct exploratory spatial data analysis. This notebook will provide an overview of ways to conduct exploratory spatial analysis in Python. 
First, let's read in some data:
End of explanation
"""


I_HR90 = ps.Moran(data.HR90.values, W)

I_HR90.I, I_HR90.p_sim

"""
Explanation: In PySAL, commonly-used analysis methods are very easy to access. For example, if we were interested in examining the spatial dependence in HR90 we could quickly compute a Moran's $I$ statistic:
End of explanation
"""


I_HR90.sim[0:5]

"""
Explanation: Thus, the $I$ statistic is $.383$ for this data, and has a very small $p$ value. 
We can visualize the distribution of simulated $I$ statistics using the stored collection of simulated statistics:
End of explanation
"""


import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

sns.kdeplot(I_HR90.sim, shade=True)
plt.vlines(I_HR90.sim, 0, 1)
plt.vlines(I_HR90.I, 0, 40, 'r')

"""
Explanation: A simple way to visualize this distribution is to make a KDEplot (like we've done before), and add a rug showing all of the simulated points, and a vertical line denoting the observed value of the statistic:
End of explanation
"""


sns.kdeplot(I_HR90.sim, shade=True)
plt.vlines(I_HR90.sim, 0, 1)
plt.vlines(I_HR90.EI+.01, 0, 40, 'r')

"""
Explanation: Instead, if our $I$ statistic were close to our expected value, I_HR90.EI, our plot might look like this:
End of explanation
"""


c_HR90 = ps.Geary(data.HR90.values, W)
#ps.Gamma
#ps.Join_Counts

c_HR90.C, c_HR90.p_sim

"""
Explanation: In addition to univariate Moran's $I$, PySAL provides many other types of autocorrelation statistics:
End of explanation
"""


bv_HRBLK = ps.Moran_BV(data.HR90.values, data.BLK90.values, W)

bv_HRBLK.I, bv_HRBLK.p_sim

"""
Explanation: Since the statistic is below one with a significant $p$-value, it indicates the same thing as the Moran's $I$ above, moderate significant global spatial dependence in HR90. 
In addition, we can compute a global Bivariate Moran statistic, which relates an observation to the spatial lag of another observation:
End of explanation
"""


LMo_HR90 = ps.Moran_Local(data.HR90.values, W)

"""
Explanation: Local Autocorrelation Statistics
In addition to the Global autocorrelation statistics, PySAL has many local autocorrelation statistics. Let's compute a local Moran statistic for the same data shown above:
End of explanation
"""


LMo_HR90.Is, LMo_HR90.p_sim

"""
Explanation: Now, instead of a single $I$ statistic, we have an array of local $I_i$ statistics, stored in the .Is attribute, and p-values from the simulation are in p_sim.
End of explanation
"""


LMo_HR90 = ps.Moran_Local(data.HR90.values, W, permutations=9999)

"""
Explanation: We can adjust the number of permutations used to derive every pseudo-$p$ value by passing a different permutations argument:
End of explanation
"""


Lag_HR90 = ps.lag_spatial(W, data.HR90.values)
HR90 = data.HR90.values

"""
Explanation: In addition to the typical clustermap, a helpful visualization for LISA statistics is a Moran scatterplot with statistically significant LISA values highlighted. 
This is very simple, if we use the same strategy we used before:
First, construct the spatial lag of the covariate:
End of explanation
"""


sigs = HR90[LMo_HR90.p_sim <= .001]
W_sigs = Lag_HR90[LMo_HR90.p_sim <= .001]
insigs = HR90[LMo_HR90.p_sim > .001]
W_insigs = Lag_HR90[LMo_HR90.p_sim > .001]

"""
Explanation: Then, we want to plot the statistically-significant LISA values in a different color than the others. To do this, first find all of the statistically significant LISAs. Since the $p$-values are in the same order as the $I_i$ statistics, we can do this in the following way
End of explanation
"""


b,a = np.polyfit(HR90, Lag_HR90, 1)

"""
Explanation: Then, since we have a lot of points, we can plot the points with a statistically insignficant LISA value lighter using the alpha keyword. In addition, we would like to plot the statistically significant points in a dark red color.
End of explanation
"""


plt.plot(sigs, W_sigs, '.', color='firebrick')
plt.plot(insigs, W_insigs, '.k', alpha=.2)
 # dashed vert at mean of the last year's PCI
plt.vlines(HR90.mean(), Lag_HR90.min(), Lag_HR90.max(), linestyle='--')
 # dashed horizontal at mean of lagged PCI
plt.hlines(Lag_HR90.mean(), HR90.min(), HR90.max(), linestyle='--')

# red line of best fit using global I as slope
plt.plot(HR90, a + b*HR90, 'r')
plt.text(s='$I = %.3f$' % I_HR90.I, x=50, y=15, fontsize=18)
plt.title('Moran Scatterplot')
plt.ylabel('Spatial Lag of HR90')
plt.xlabel('HR90')

"""
Explanation: Matplotlib has a list of named colors and will interpret colors that are provided in hexadecimal strings:
End of explanation
"""


pd.get_dummies(data.SOUTH) #dummies for south (already binary)

"""
Explanation: We can also make a LISA map of the data. 
Simple exploratory regression
Sometimes, to check for simple spatial heterogeneity, a fixed effects model can be estimated.  If the heterogeneity has known bounds. First, though, note that pandas can build dummy variable matrices from categorical data very quickly:
End of explanation
"""


pd.get_dummies(data.STATE_NAME) #dummies for state by name

"""
Explanation: Where this becomes handy is if you have data you know can be turned into a dummy variable, but is not yet correctly encoded. 
For example, the same call as above can make a dummy variable matrix to encode state fixed effects using the STATE_NAME variable:
End of explanation
"""


y = data[['HR90']].values
x = data[['BLK90']].values

unrelated_effect = np.random.normal(0,100, size=y.shape[0]).reshape(y.shape)

X = np.hstack((x, unrelated_effect))

regimes = data.SOUTH.values.tolist()

regime_reg = ps.spreg.OLS_Regimes(y, X, regimes)

betas = regime_reg.betas
sebetas = np.sqrt(regime_reg.vm.diagonal())

sebetas

plt.plot(betas,'ok')
plt.axis([-1,6,-1,8])
plt.hlines(0,-1,6, color='k', linestyle='--')
plt.errorbar([0,1,2,3,4,5], betas.flatten(), yerr=sebetas*3, fmt='o', ecolor='r')
plt.xticks([-.5,0.5,1.5,2.5,3.5,4.5, 5.5], ['',
                                       'Not South: Constant',
                                       'Not South: BLK90',
                                       'Not South: Not South',
                                       'South: Constant',
                                       'South: South',
                                       'South: Unrelated',''], rotation='325')
plt.title('Regime Fixed Effects')

"""
Explanation: For now, let's estimate a spatial regimes regression on the south/not-south division. To show how a regimes effects plot may look, let's consider one covariate that is likely related and one that is very likely unrelated to $y$. That is our formal specification for the regression will be:
$$ y = \beta_{0} + x_{1a}\beta_{1a} + x_{1b}\beta_{1b} + x_{2}\beta_{2} + \epsilon $$
Where $x_{1a} = 1$ when an observation is not in the south, and zero when it is not. This is a simple spatial fixed effects setup, where each different spatial unit is treated as having a special effect on observations inside of it. 
In addition, we'll add an unrelated term to show how the regression effects visualization works:
End of explanation
"""


