In [920]:
import numpy as np
import pandas as pd
import os
os.chdir('/Users/toshan/dev/pysal/pysal/contrib/spint')
from gravity import Gravity, Production, Attraction, Doubly, BaseGravity
import statsmodels.formula.api as smf
from statsmodels.api import families
import pysal.spreg.ml_lag as ml_lag
import pysal
os.chdir('/Users/toshan/dev/pysal/pysal/weights')
from spintW import ODW
from pysal.spreg.twosls_sp import GM_Lag, BaseGM_Lag
from pysal.spreg.twosls import BaseTSLS
from pysal.spreg.utils import set_endog
In [921]:
def m_inverse_prod(w, data, scalar, post_multiply=False, inv_method="power_exp", threshold=0.0000000001, max_iterations=None):
"""
Parameters
----------
w : Pysal W object
nxn Pysal spatial weights object
data : Numpy array
nx1 vector of data
scalar : float
Scalar value (typically rho or lambda)
post_multiply : boolean
If True then post-multiplies the data vector by the
inverse of the spatial filter, if false then
pre-multiplies.
threshold : float
Test value to stop the iterations. Test is against
sqrt(increment' * increment), where increment is a
vector representing the contribution from each
iteration.
max_iterations : integer
Maximum number of iterations for the expansion.
Examples
--------
"""
N = data.shape[0]
matrix = la.inv(np.eye(N))
try:
for i in range(len(w)):
matrix -= (scalar[i] * w[i].full()[0])
except:
for i in range(len(w)):
matrix -= (scalar[i] * w[i])
if post_multiply:
inv_prod = spdot(data.T, np.asarray(matrix))
else:
inv_prod = spdot(np.asarray(matrix), data)
return inv_prod
import numpy as np
import numpy.linalg as la
from scipy import sparse as sp
from scipy.sparse.linalg import splu as SuperLU
import pysal as ps
from pysal.spreg.utils import RegressionPropsY, RegressionPropsVM, inverse_prod, spdot
import pysal.spreg.diagnostics as DIAG
import pysal.spreg.user_output as USER
import pysal.spreg.summary_output as SUMMARY
from pysal.spreg.w_utils import symmetrize
try:
from scipy.optimize import minimize_scalar
from scipy.optimize import minimize
minimize_scalar_available = True
except ImportError:
minimize_scalar_available = False
__all__ = ["ML_Lag"]
class BaseML_Lag(RegressionPropsY, RegressionPropsVM):
"""
ML estimation of the spatial lag model (note no consistency
checks, diagnostics or constants added); Anselin (1988) [Anselin1988]_
Parameters
----------
y : array
nx1 array for dependent variable
x : array
Two dimensional array with n rows and one column for each
independent (exogenous) variable, excluding the constant
w : pysal W object
Spatial weights object
method : string
if 'full', brute force calculation (full matrix expressions)
if 'ord', Ord eigenvalue method
if 'LU', LU sparse matrix decomposition
epsilon : float
tolerance criterion in mimimize_scalar function and inverse_product
Attributes
----------
betas : array
(k+1)x1 array of estimated coefficients (rho first)
rho : float
estimate of spatial autoregressive coefficient
u : array
nx1 array of residuals
predy : array
nx1 array of predicted y values
n : integer
Number of observations
k : integer
Number of variables for which coefficients are estimated
(including the constant, excluding the rho)
y : array
nx1 array for dependent variable
x : array
Two dimensional array with n rows and one column for each
independent (exogenous) variable, including the constant
method : string
log Jacobian method
if 'full': brute force (full matrix computations)
if 'ord' : Ord eigenvalue method
epsilon : float
tolerance criterion used in minimize_scalar function and inverse_product
mean_y : float
Mean of dependent variable
std_y : float
Standard deviation of dependent variable
vm : array
Variance covariance matrix (k+1 x k+1)
vm1 : array
Variance covariance matrix (k+2 x k+2) includes sigma2
sig2 : float
Sigma squared used in computations
logll : float
maximized log-likelihood (including constant terms)
predy_e : array
predicted values from reduced form
e_pred : array
prediction errors using reduced form predicted values
Examples
--------
>>> import numpy as np
>>> import pysal as ps
>>> db = ps.open(ps.examples.get_path("baltim.dbf"),'r')
>>> ds_name = "baltim.dbf"
>>> y_name = "PRICE"
>>> y = np.array(db.by_col(y_name)).T
>>> y.shape = (len(y),1)
>>> x_names = ["NROOM","NBATH","PATIO","FIREPL","AC","GAR","AGE","LOTSZ","SQFT"]
>>> x = np.array([db.by_col(var) for var in x_names]).T
>>> x = np.hstack((np.ones((len(y),1)),x))
>>> ww = ps.open(ps.examples.get_path("baltim_q.gal"))
>>> w = ww.read()
>>> ww.close()
>>> w.transform = 'r'
>>> w_name = "baltim_q.gal"
>>> mllag = BaseML_Lag(y,x,w,method='ord') #doctest: +SKIP
>>> "{0:.6f}".format(mllag.rho) #doctest: +SKIP
'0.425885'
>>> np.around(mllag.betas, decimals=4) #doctest: +SKIP
array([[ 4.3675],
[ 0.7502],
[ 5.6116],
[ 7.0497],
[ 7.7246],
[ 6.1231],
[ 4.6375],
[-0.1107],
[ 0.0679],
[ 0.0794],
[ 0.4259]])
>>> "{0:.6f}".format(mllag.mean_y) #doctest: +SKIP
'44.307180'
>>> "{0:.6f}".format(mllag.std_y) #doctest: +SKIP
'23.606077'
>>> np.around(np.diag(mllag.vm1), decimals=4) #doctest: +SKIP
array([ 23.8716, 1.1222, 3.0593, 7.3416, 5.6695, 5.4698,
2.8684, 0.0026, 0.0002, 0.0266, 0.0032, 220.1292])
>>> np.around(np.diag(mllag.vm), decimals=4) #doctest: +SKIP
array([ 23.8716, 1.1222, 3.0593, 7.3416, 5.6695, 5.4698,
2.8684, 0.0026, 0.0002, 0.0266, 0.0032])
>>> "{0:.6f}".format(mllag.sig2) #doctest: +SKIP
'151.458698'
>>> "{0:.6f}".format(mllag.logll) #doctest: +SKIP
'-832.937174'
>>> mllag = BaseML_Lag(y,x,w) #doctest: +SKIP
>>> "{0:.6f}".format(mllag.rho) #doctest: +SKIP
'0.425885'
>>> np.around(mllag.betas, decimals=4) #doctest: +SKIP
array([[ 4.3675],
[ 0.7502],
[ 5.6116],
[ 7.0497],
[ 7.7246],
[ 6.1231],
[ 4.6375],
[-0.1107],
[ 0.0679],
[ 0.0794],
[ 0.4259]])
>>> "{0:.6f}".format(mllag.mean_y) #doctest: +SKIP
'44.307180'
>>> "{0:.6f}".format(mllag.std_y) #doctest: +SKIP
'23.606077'
>>> np.around(np.diag(mllag.vm1), decimals=4) #doctest: +SKIP
array([ 23.8716, 1.1222, 3.0593, 7.3416, 5.6695, 5.4698,
2.8684, 0.0026, 0.0002, 0.0266, 0.0032, 220.1292])
>>> np.around(np.diag(mllag.vm), decimals=4) #doctest: +SKIP
array([ 23.8716, 1.1222, 3.0593, 7.3416, 5.6695, 5.4698,
2.8684, 0.0026, 0.0002, 0.0266, 0.0032])
>>> "{0:.6f}".format(mllag.sig2) #doctest: +SKIP
'151.458698'
>>> "{0:.6f}".format(mllag.logll) #doctest: +SKIP
'-832.937174'
"""
def __init__(self, y, x, wo, wd, ww, method='full', epsilon=0.0000001):
# set up main regression variables and spatial filters
self.y = y
self.x = x
self.n, self.k = self.x.shape
self.method = method
self.epsilon = epsilon
#W = w.full()[0]
#Wsp = w.sparse
oylag = ps.lag_spatial(wo, y)
dylag = ps.lag_spatial(wd, y)
wylag = ps.lag_spatial(ww, y)
ylag = (wo.sparse*wd.sparse*ODw.sparse)*f
# b0, b1, e0 and e1
xtx = spdot(self.x.T, self.x)
xtxi = la.inv(xtx)
xty = spdot(self.x.T, self.y)
xtylo = spdot(self.x.T, oylag)
xtyld = spdot(self.x.T, dylag)
xtylw = spdot(self.x.T, wylag)
xtyl = spdot(self.x.T, ylag)
b0 = np.dot(xtxi, xty)
b1 = np.dot(xtxi, xtylo)
b2 = np.dot(xtxi, xtyld)
b3 = np.dot(xtxi, xtylw)
b4 = np.dot(xtxi, xtyl)
e0 = self.y - spdot(x, b0)
e1 = oylag - spdot(x, b1)
e2 = dylag - spdot(x, b2)
e3 = wylag - spdot(x, b3)
e4 = ylag - spdot(x, b4)
methodML = method.upper()
# call minimizer using concentrated log-likelihood to get rho
if methodML in ['FULL', 'LU', 'ORD']:
if methodML == 'FULL':
Wo = wo.full()[0] # moved here
Wd = wd.full()[0]
Ww = ww.full()[0]
x0 = np.array([0.0, 0.0, 0.0])
bounds = [(-1.0, 1.0), (-1.0, 1.0), (-1.0, 1.0)]
res = minimize(lag_c_loglik, x0 , bounds=bounds,
args=(
self.n, e0, e1, e2, e3, e4, Wo, Wd, Ww),
tol=epsilon)
elif methodML == 'LU':
I = sp.identity(w.n)
Wsp = w.sparse # moved here
res = minimize_scalar(lag_c_loglik_sp, 0.0, bounds=(-1.0,1.0),
args=(self.n, e0, e1, I, Wsp),
method='bounded', tol=epsilon)
elif methodML == 'ORD':
# check on symmetry structure
if w.asymmetry(intrinsic=False) == []:
ww = symmetrize(w)
WW = ww.todense()
evals = la.eigvalsh(WW)
else:
W = w.full()[0] # moved here
evals = la.eigvals(W)
res = minimize_scalar(lag_c_loglik_ord, 0.0, bounds=(-1.0, 1.0),
args=(
self.n, e0, e1, evals), method='bounded',
tol=epsilon)
else:
# program will crash, need to catch
print("{0} is an unsupported method".format(methodML))
self = None
return
#print res.x
self.o_rho, self.d_rho, self.w_rho = res.x
# compute full log-likelihood, including constants
ln2pi = np.log(2.0 * np.pi)
llik = -res.fun - self.n / 2.0 * ln2pi - self.n / 2.0
self.logll = llik[0][0]
# b, residuals and predicted values
b = b0 - self.o_rho * b1 - self.d_rho * b2 - self.w_rho *b3
self.betas = np.vstack((b, self.o_rho, self.d_rho, self.w_rho)) # rho added as last coefficient
self.u = e0 - self.o_rho * e1 - self.o_rho * e2 - self.w_rho * e3
self.predy = self.y - self.u
xb = spdot(x, b)
self.predy_e = m_inverse_prod(
[wo.sparse, wd.sparse, ww.sparse], xb, [self.o_rho, self.d_rho, self.w_rho], inv_method="power_exp", threshold=epsilon)
self.e_pred = self.y - self.predy_e
# residual variance
self.sig2 = self.sig2n # no allowance for division by n-k
#print self.betas
# information matrix
a = -self.o_rho * Wo - self.d_rho * Wd + self.w_rho * Ww
W_o = -Wo - self.d_rho * Wd + self.w_rho * Ww
W_d = -self.o_rho * Wo - Wd + self.w_rho * Ww
W_w = -self.o_rho * Wo - self.d_rho * Wd + Ww
np.fill_diagonal(a, 1.0)
ai = la.inv(a)
o_wai = np.dot(W_o, ai)
o_tr1 = np.trace(o_wai)
d_wai = np.dot(W_d, ai)
d_tr1 = np.trace(d_wai)
w_wai = np.dot(W_w, ai)
w_tr1 = np.trace(w_wai)
oo_wai2 = np.dot(o_wai, o_wai)
oo_tr2 = np.trace(oo_wai2)
dd_wai2 = np.dot(d_wai, d_wai)
dd_tr2 = np.trace(dd_wai2)
ww_wai2 = np.dot(w_wai, w_wai)
ww_tr2 = np.trace(ww_wai2)
od_wai2 = np.dot(o_wai, d_wai)
od_tr2 = np.trace(od_wai2)
do_wai2 = np.dot(d_wai, o_wai)
do_tr2 = np.trace(do_wai2)
ow_wai2 = np.dot(o_wai, w_wai)
ow_tr2 = np.trace(ow_wai2)
wo_wai2 = np.dot(w_wai, o_wai)
wo_tr2 = np.trace(wo_wai2)
dw_wai2 = np.dot(d_wai, w_wai)
dw_tr2 = np.trace(dw_wai2)
wd_wai2 = np.dot(w_wai, d_wai)
wd_tr2 = np.trace(wd_wai2)
oo_waiTwai = np.dot(o_wai.T, o_wai)
oo_tr3 = np.trace(oo_waiTwai)
dd_waiTwai = np.dot(d_wai.T, d_wai)
dd_tr3 = np.trace(dd_waiTwai)
ww_waiTwai = np.dot(w_wai.T, w_wai)
ww_tr3 = np.trace(ww_waiTwai)
od_waiTwai = np.dot(o_wai.T, d_wai)
od_tr3 = np.trace(od_waiTwai)
do_waiTwai = np.dot(d_wai.T, o_wai)
do_tr3 = np.trace(do_waiTwai)
ow_waiTwai = np.dot(o_wai.T, w_wai)
ow_tr3 = np.trace(ow_waiTwai)
wo_waiTwai = np.dot(w_wai.T, o_wai)
wo_tr3 = np.trace(wo_waiTwai)
dw_waiTwai = np.dot(d_wai.T, w_wai)
dw_tr3 = np.trace(dw_waiTwai)
wd_waiTwai = np.dot(w_wai.T, d_wai)
wd_tr3 = np.trace(wd_waiTwai)
o_wpredy = ps.lag_spatial(wo, self.predy_e)
d_wpredy = ps.lag_spatial(wd, self.predy_e)
w_wpredy = ps.lag_spatial(ww, self.predy_e)
o_xTwpy = spdot(x.T, o_wpredy)
d_xTwpy = spdot(x.T, d_wpredy)
w_xTwpy = spdot(x.T, w_wpredy)
oo_wpyTwpy = np.dot(o_wpredy.T, o_wpredy)
dd_wpyTwpy = np.dot(d_wpredy.T, d_wpredy)
ww_wpyTwpy = np.dot(w_wpredy.T, w_wpredy)
od_wpyTwpy = np.dot(o_wpredy.T, d_wpredy)
do_wpyTwpy = np.dot(d_wpredy.T, o_wpredy)
ow_wpyTwpy = np.dot(o_wpredy.T, w_wpredy)
wo_wpyTwpy = np.dot(w_wpredy.T, o_wpredy)
dw_wpyTwpy = np.dot(d_wpredy.T, w_wpredy)
wd_wpyTwpy = np.dot(w_wpredy.T, d_wpredy)
sig2 = self.sig2
beta_beta = xtx / sig2
beta_o_rho = o_xTwpy.T / sig2
beta_d_rho = d_xTwpy.T / sig2
beta_w_rho = w_xTwpy.T / sig2
beta_sig2 = np.zeros((1, self.k))
o_rho_beta = o_xTwpy / sig2
o_rho_o_rho = oo_tr2 + oo_tr3 + oo_wpyTwpy / sig2
o_rho_d_rho = od_tr2 + od_tr3 + od_wpyTwpy / sig2
o_rho_w_rho = ow_tr2 + ow_tr3 + ow_wpyTwpy / sig2
o_rho_sig2 = o_tr1 / sig2
d_rho_beta = d_xTwpy / sig2
d_rho_o_rho = do_tr2 + do_tr3 + do_wpyTwpy / sig2
d_rho_d_rho = dd_tr2 + dd_tr3 + dd_wpyTwpy / sig2
d_rho_w_rho = dw_tr2 + dw_tr3 + dw_wpyTwpy / sig2
d_rho_sig2 = d_tr1 / sig2
w_rho_beta = w_xTwpy / sig2
w_rho_o_rho = wo_tr2 + wo_tr3 + wo_wpyTwpy / sig2
w_rho_d_rho = wd_tr2 + wd_tr3 + wd_wpyTwpy / sig2
w_rho_w_rho = ww_tr2 + ww_tr3 + ww_wpyTwpy / sig2
w_rho_sig2 = w_tr1 / sig2
sig2_beta = np.zeros((self.k, 1))
sig2_o_rho = o_tr1 / sig2
sig2_d_rho = d_tr1 / sig2
sig2_w_rho = w_tr1 / sig2
sig2_sig2 = self.n / (2.0 * sig2 ** 2)
# order of variables is beta, rho_o, rho_d, rho_w, sigma2
v1 = np.vstack((beta_beta, beta_o_rho, beta_d_rho, beta_w_rho, beta_sig2))
v2 = np.vstack((o_rho_beta, o_rho_o_rho, o_rho_d_rho, o_rho_w_rho, o_rho_sig2))
v3 = np.vstack((d_rho_beta, d_rho_o_rho, d_rho_d_rho, d_rho_w_rho, d_rho_sig2))
v4 = np.vstack((w_rho_beta, w_rho_d_rho, w_rho_d_rho, w_rho_w_rho, w_rho_sig2))
v5 = np.vstack((sig2_beta, sig2_d_rho, sig2_d_rho, sig2_w_rho, sig2_sig2))
v = np.hstack((v1, v2, v3, v4, v5))
self.vm1 = la.inv(v) # vm1 includes variance for sigma2
self.vm = self.vm1[:-1, :-1] # vm is for coefficients only
def lag_c_loglik(rho, n, e0, e1, e2, e3, e4, Wo, Wd, Ww):
# concentrated log-lik for lag model, no constants, brute force
o_rho, d_rho, w_rho = rho
er = e0 - (o_rho * e1) - (d_rho * e2) - (w_rho * e3)
sig2 = np.dot(er.T, er) / n
nlsig2 = (n / 2.0) * np.log(sig2)
a = (np.identity(n) - ((o_rho * Wo)*(d_rho * Wd)*(w_rho * Ww)))
np.fill_diagonal(a, 1.0)
jacob = np.log(np.linalg.det(a))
# this is the negative of the concentrated log lik for minimization
clik = nlsig2 - jacob
#print clik
return clik
def lag_c_loglik_sp(rho, n, e0, e1, I, Wsp):
# concentrated log-lik for lag model, sparse algebra
if isinstance(rho, np.ndarray):
if rho.shape == (1,1):
rho = rho[0][0] #why does the interior value change?
er = e0 - rho * e1
sig2 = np.dot(er.T, er) / n
nlsig2 = (n / 2.0) * np.log(sig2)
a = I - rho * Wsp
LU = SuperLU(a.tocsc())
jacob = np.sum(np.log(np.abs(LU.U.diagonal())))
clike = nlsig2 - jacob
return clike
def lag_c_loglik_ord(rho, n, e0, e1, evals):
# concentrated log-lik for lag model, no constants, Ord eigenvalue method
er = e0 - rho * e1
sig2 = np.dot(er.T, er) / n
nlsig2 = (n / 2.0) * np.log(sig2)
revals = rho * evals
jacob = np.log(1 - revals).sum()
if isinstance(jacob, complex):
jacob = jacob.real
# this is the negative of the concentrated log lik for minimization
clik = nlsig2 - jacob
return clik
SAR SI DGP
$y = (I - \rho W) X \beta + (I - \rho W)\epsilon$
$X \beta = X_{o} \beta_{o} + X_{d} \beta_{d} + X_{d_{ij}} \beta_{d_{ij}}$
$\rho W = (I - \rho_{o}W_{o})(I - \rho_{d}W_{d}) = (I - \rho_{o}W_{o} - \rho_{d}W_{d} - \rho_{w}W_{w})$
$\rho_{w}W_{w} = \rho_{o} \rho_{d}W_{o}*W_{d}$
SAR SI Model Specification
$y = \rho_{o}W_{o} + \rho_{d}W_{d} + \rho_{w}W_{w} + k + X \beta + (I - \rho W)\epsilon$
In [922]:
#Simulate covariates
o_vars = np.random.randint(5000,10000, (81,1))
d_vars = np.random.randint(5000,10000, (81,1))
dij = np.random.randint(50,1000, (81,1))
X = np.hstack([o_vars.reshape((-1,1)), d_vars.reshape((-1,1)), dij.reshape((-1,1))])
In [954]:
#Simulate spatial weights and compute dependent variable using DGP
z = 7
n = 49
N = n**2
o_vars = np.random.randint(5000,10000, (N,1))
d_vars = np.random.randint(5000,10000, (N,1))
dij = np.random.randint(50,1000, (N,1))
X = np.hstack([o_vars.reshape((-1,1)), d_vars.reshape((-1,1)), dij.reshape((-1,1))])
o = pysal.weights.lat2W(z, z)
d = pysal.weights.lat2W(z, z)
wo = np.kron(o.full()[0], np.identity(n))
wo = pysal.weights.full2W(wo)
wo.transform = 'r'
wd = np.kron(np.identity(n), d.full()[0])
wd = pysal.weights.full2W(wd)
wd.transform = 'r'
ODw = ODW(o,d)
ODw.transform = 'r'
o_rho = .15
d_rho = .15
w_rho = -.3
ar = np.array(np.linalg.inv((np.identity(N) - o_rho*wo.sparse - d_rho*wd.sparse - w_rho*ODw.sparse)))
betas = np.vstack([1,.3,.3,-5])
X_constant = np.hstack([np.ones((N,1)), X])
XB = np.dot(X_constant,betas)
e = np.random.normal(0, 100, (N,1))
y = np.dot(ar, XB) + np.dot(ar, e)
In [955]:
#Prep instruments
w_lags=1
o_yend, o_q = set_endog(y, X, wo, None, None, w_lags, True)
d_yend, d_q = set_endog(y, X, wd, None, None, w_lags, True)
w_yend, w_q = set_endog(y, X, ODw, None, None, w_lags, True)
In [959]:
#GM_Lag via TSLS
w = (wo.sparse, wd.sparse, ODw.sparse)
yend = np.hstack([o_yend, d_yend, w_yend])
q = np.hstack([o_q, d_q, w_q])
model = BaseTSLS(y=y, x=X_constant, yend=yend, q=q)
print model.betas
print model.vm.diagonal()
In [958]:
#Extended ML_Lag
model = BaseML_Lag(y, X_constant, wo, wd, ODw)
print model.betas
print model.vm.diagonal()
In [964]:
print len(y>=0)
In [965]:
len(y)
Out[965]:
In [ ]:
generate 2000