1""" 2ML Estimation of Spatial Lag Model 3""" 4 5__author__ = "Luc Anselin luc.anselin@asu.edu, \ 6 Serge Rey srey@asu.edu, \ 7 Levi Wolf levi.john.wolf@gmail.com" 8 9import numpy as np 10import numpy.linalg as la 11from scipy import sparse as sp 12from scipy.sparse.linalg import splu as SuperLU 13from .utils import RegressionPropsY, RegressionPropsVM, inverse_prod, set_warn 14from .sputils import spdot, spfill_diagonal, spinv, spbroadcast 15from . import diagnostics as DIAG 16from . import user_output as USER 17from . import summary_output as SUMMARY 18from .w_utils import symmetrize 19from libpysal import weights 20try: 21 from scipy.optimize import minimize_scalar 22 minimize_scalar_available = True 23except ImportError: 24 minimize_scalar_available = False 25 26__all__ = ["ML_Lag"] 27 28 29class BaseML_Lag(RegressionPropsY, RegressionPropsVM): 30 31 """ 32 ML estimation of the spatial lag model (note no consistency 33 checks, diagnostics or constants added) :cite:`Anselin1988` 34 35 Parameters 36 ---------- 37 y : array 38 nx1 array for dependent variable 39 x : array 40 Two dimensional array with n rows and one column for each 41 independent (exogenous) variable, excluding the constant 42 w : pysal W object 43 Spatial weights object 44 method : string 45 if 'full', brute force calculation (full matrix expressions) 46 if 'ord', Ord eigenvalue method 47 if 'LU', LU sparse matrix decomposition 48 epsilon : float 49 tolerance criterion in mimimize_scalar function and inverse_product 50 51 Attributes 52 ---------- 53 betas : array 54 (k+1)x1 array of estimated coefficients (rho first) 55 rho : float 56 estimate of spatial autoregressive coefficient 57 u : array 58 nx1 array of residuals 59 predy : array 60 nx1 array of predicted y values 61 n : integer 62 Number of observations 63 k : integer 64 Number of variables for which coefficients are estimated 65 (including the constant, excluding the rho) 66 y : array 67 nx1 array for dependent variable 68 x : array 69 Two dimensional array with n rows and one column for each 70 independent (exogenous) variable, including the constant 71 method : string 72 log Jacobian method 73 if 'full': brute force (full matrix computations) 74 if 'ord' : Ord eigenvalue method 75 epsilon : float 76 tolerance criterion used in minimize_scalar function and inverse_product 77 mean_y : float 78 Mean of dependent variable 79 std_y : float 80 Standard deviation of dependent variable 81 vm : array 82 Variance covariance matrix (k+1 x k+1) 83 vm1 : array 84 Variance covariance matrix (k+2 x k+2) includes sigma2 85 sig2 : float 86 Sigma squared used in computations 87 logll : float 88 maximized log-likelihood (including constant terms) 89 predy_e : array 90 predicted values from reduced form 91 e_pred : array 92 prediction errors using reduced form predicted values 93 94 95 Examples 96 -------- 97 98 >>> import numpy as np 99 >>> import libpysal 100 >>> from libpysal.examples import load_example 101 >>> import geopandas as gpd 102 >>> from libpysal.weights import Queen 103 >>> import spreg 104 >>> np.set_printoptions(suppress=True) #prevent scientific format 105 >>> baltimore = load_example('Baltimore') 106 >>> db = libpysal.io.open(baltimore.get_path("baltim.dbf"),'r') 107 >>> df = gpd.read_file(baltimore.get_path("baltim.shp")) 108 >>> ds_name = "baltim.dbf" 109 >>> y_name = "PRICE" 110 >>> y = np.array(db.by_col(y_name)).T 111 >>> y.shape = (len(y),1) 112 >>> x_names = ["NROOM","NBATH","PATIO","FIREPL","AC","GAR","AGE","LOTSZ","SQFT"] 113 >>> x = np.array([db.by_col(var) for var in x_names]).T 114 >>> x = np.hstack((np.ones((len(y),1)),x)) 115 >>> w = Queen.from_dataframe(df) 116 >>> w.transform = 'r' 117 >>> w_name = "baltim_q.gal" 118 >>> mllag = spreg.ml_lag.BaseML_Lag(y,x,w,method='ord') #doctest: +SKIP 119 >>> "{0:.6f}".format(mllag.rho) #doctest: +SKIP 120 '0.425885' 121 >>> np.around(mllag.betas, decimals=4) #doctest: +SKIP 122 array([[ 4.3675], 123 [ 0.7502], 124 [ 5.6116], 125 [ 7.0497], 126 [ 7.7246], 127 [ 6.1231], 128 [ 4.6375], 129 [-0.1107], 130 [ 0.0679], 131 [ 0.0794], 132 [ 0.4259]]) 133 >>> "{0:.6f}".format(mllag.mean_y) #doctest: +SKIP 134 '44.307180' 135 >>> "{0:.6f}".format(mllag.std_y) #doctest: +SKIP 136 '23.606077' 137 >>> np.around(np.diag(mllag.vm1), decimals=4) #doctest: +SKIP 138 array([ 23.8716, 1.1222, 3.0593, 7.3416, 5.6695, 5.4698, 139 2.8684, 0.0026, 0.0002, 0.0266, 0.0032, 220.1292]) 140 >>> np.around(np.diag(mllag.vm), decimals=4) #doctest: +SKIP 141 array([ 23.8716, 1.1222, 3.0593, 7.3416, 5.6695, 5.4698, 142 2.8684, 0.0026, 0.0002, 0.0266, 0.0032]) 143 >>> "{0:.6f}".format(mllag.sig2) #doctest: +SKIP 144 '151.458698' 145 >>> "{0:.6f}".format(mllag.logll) #doctest: +SKIP 146 '-832.937174' 147 >>> mllag = spreg.ml_lag.BaseML_Lag(y,x,w) #doctest: +SKIP 148 >>> "{0:.6f}".format(mllag.rho) #doctest: +SKIP 149 '0.425885' 150 >>> np.around(mllag.betas, decimals=4) #doctest: +SKIP 151 array([[ 4.3675], 152 [ 0.7502], 153 [ 5.6116], 154 [ 7.0497], 155 [ 7.7246], 156 [ 6.1231], 157 [ 4.6375], 158 [-0.1107], 159 [ 0.0679], 160 [ 0.0794], 161 [ 0.4259]]) 162 >>> "{0:.6f}".format(mllag.mean_y) #doctest: +SKIP 163 '44.307180' 164 >>> "{0:.6f}".format(mllag.std_y) #doctest: +SKIP 165 '23.606077' 166 >>> np.around(np.diag(mllag.vm1), decimals=4) #doctest: +SKIP 167 array([ 23.8716, 1.1222, 3.0593, 7.3416, 5.6695, 5.4698, 168 2.8684, 0.0026, 0.0002, 0.0266, 0.0032, 220.1292]) 169 >>> np.around(np.diag(mllag.vm), decimals=4) #doctest: +SKIP 170 array([ 23.8716, 1.1222, 3.0593, 7.3416, 5.6695, 5.4698, 171 2.8684, 0.0026, 0.0002, 0.0266, 0.0032]) 172 >>> "{0:.6f}".format(mllag.sig2) #doctest: +SKIP 173 '151.458698' 174 >>> "{0:.6f}".format(mllag.logll) #doctest: +SKIP 175 '-832.937174' 176 177 178 """ 179 180 def __init__(self, y, x, w, method='full', epsilon=0.0000001): 181 # set up main regression variables and spatial filters 182 self.y = y 183 self.x = x 184 self.n, self.k = self.x.shape 185 self.method = method 186 self.epsilon = epsilon 187 #W = w.full()[0] 188 #Wsp = w.sparse 189 ylag = weights.lag_spatial(w, y) 190 # b0, b1, e0 and e1 191 xtx = spdot(self.x.T, self.x) 192 xtxi = la.inv(xtx) 193 xty = spdot(self.x.T, self.y) 194 xtyl = spdot(self.x.T, ylag) 195 b0 = spdot(xtxi, xty) 196 b1 = spdot(xtxi, xtyl) 197 e0 = self.y - spdot(x, b0) 198 e1 = ylag - spdot(x, b1) 199 methodML = method.upper() 200 # call minimizer using concentrated log-likelihood to get rho 201 if methodML in ['FULL', 'LU', 'ORD']: 202 if methodML == 'FULL': 203 W = w.full()[0] # moved here 204 res = minimize_scalar(lag_c_loglik, 0.0, bounds=(-1.0, 1.0), 205 args=( 206 self.n, e0, e1, W), method='bounded', 207 tol=epsilon) 208 elif methodML == 'LU': 209 I = sp.identity(w.n) 210 Wsp = w.sparse # moved here 211 W = Wsp 212 res = minimize_scalar(lag_c_loglik_sp, 0.0, bounds=(-1.0,1.0), 213 args=(self.n, e0, e1, I, Wsp), 214 method='bounded', tol=epsilon) 215 elif methodML == 'ORD': 216 # check on symmetry structure 217 if w.asymmetry(intrinsic=False) == []: 218 ww = symmetrize(w) 219 WW = np.array(ww.todense()) 220 evals = la.eigvalsh(WW) 221 W = WW 222 else: 223 W = w.full()[0] # moved here 224 evals = la.eigvals(W) 225 res = minimize_scalar(lag_c_loglik_ord, 0.0, bounds=(-1.0, 1.0), 226 args=( 227 self.n, e0, e1, evals), method='bounded', 228 tol=epsilon) 229 else: 230 # program will crash, need to catch 231 print(("{0} is an unsupported method".format(methodML))) 232 self = None 233 return 234 235 self.rho = res.x[0][0] 236 237 # compute full log-likelihood, including constants 238 ln2pi = np.log(2.0 * np.pi) 239 llik = -res.fun - self.n / 2.0 * ln2pi - self.n / 2.0 240 self.logll = llik[0][0] 241 242 # b, residuals and predicted values 243 244 b = b0 - self.rho * b1 245 self.betas = np.vstack((b, self.rho)) # rho added as last coefficient 246 self.u = e0 - self.rho * e1 247 self.predy = self.y - self.u 248 249 xb = spdot(x, b) 250 251 self.predy_e = inverse_prod( 252 w.sparse, xb, self.rho, inv_method="power_exp", threshold=epsilon) 253 self.e_pred = self.y - self.predy_e 254 255 # residual variance 256 self._cache = {} 257 self.sig2 = self.sig2n # no allowance for division by n-k 258 259 # information matrix 260 # if w should be kept sparse, how can we do the following: 261 a = -self.rho * W 262 spfill_diagonal(a, 1.0) 263 ai = spinv(a) 264 wai = spdot(W, ai) 265 tr1 = wai.diagonal().sum() #same for sparse and dense 266 267 wai2 = spdot(wai, wai) 268 tr2 = wai2.diagonal().sum() 269 270 waiTwai = spdot(wai.T, wai) 271 tr3 = waiTwai.diagonal().sum() 272 ### to here 273 274 wpredy = weights.lag_spatial(w, self.predy_e) 275 wpyTwpy = spdot(wpredy.T, wpredy) 276 xTwpy = spdot(x.T, wpredy) 277 278 # order of variables is beta, rho, sigma2 279 280 v1 = np.vstack( 281 (xtx / self.sig2, xTwpy.T / self.sig2, np.zeros((1, self.k)))) 282 v2 = np.vstack( 283 (xTwpy / self.sig2, tr2 + tr3 + wpyTwpy / self.sig2, tr1 / self.sig2)) 284 v3 = np.vstack( 285 (np.zeros((self.k, 1)), tr1 / self.sig2, self.n / (2.0 * self.sig2 ** 2))) 286 287 v = np.hstack((v1, v2, v3)) 288 289 self.vm1 = la.inv(v) # vm1 includes variance for sigma2 290 self.vm = self.vm1[:-1, :-1] # vm is for coefficients only 291 292 293class ML_Lag(BaseML_Lag): 294 295 """ 296 ML estimation of the spatial lag model with all results and diagnostics; :cite:`Anselin1988` 297 298 Parameters 299 ---------- 300 y : array 301 nx1 array for dependent variable 302 x : array 303 Two dimensional array with n rows and one column for each 304 independent (exogenous) variable, excluding the constant 305 w : pysal W object 306 Spatial weights object 307 method : string 308 if 'full', brute force calculation (full matrix expressions) 309 if 'ord', Ord eigenvalue method 310 epsilon : float 311 tolerance criterion in mimimize_scalar function and inverse_product 312 vm : boolean 313 if True, include variance-covariance matrix in summary 314 results 315 name_y : string 316 Name of dependent variable for use in output 317 name_x : list of strings 318 Names of independent variables for use in output 319 name_w : string 320 Name of weights matrix for use in output 321 name_ds : string 322 Name of dataset for use in output 323 324 Attributes 325 ---------- 326 betas : array 327 (k+1)x1 array of estimated coefficients (rho first) 328 rho : float 329 estimate of spatial autoregressive coefficient 330 u : array 331 nx1 array of residuals 332 predy : array 333 nx1 array of predicted y values 334 n : integer 335 Number of observations 336 k : integer 337 Number of variables for which coefficients are estimated 338 (including the constant, excluding the rho) 339 y : array 340 nx1 array for dependent variable 341 x : array 342 Two dimensional array with n rows and one column for each 343 independent (exogenous) variable, including the constant 344 method : string 345 log Jacobian method 346 if 'full': brute force (full matrix computations) 347 epsilon : float 348 tolerance criterion used in minimize_scalar function and inverse_product 349 mean_y : float 350 Mean of dependent variable 351 std_y : float 352 Standard deviation of dependent variable 353 vm : array 354 Variance covariance matrix (k+1 x k+1), all coefficients 355 vm1 : array 356 Variance covariance matrix (k+2 x k+2), includes sig2 357 sig2 : float 358 Sigma squared used in computations 359 logll : float 360 maximized log-likelihood (including constant terms) 361 aic : float 362 Akaike information criterion 363 schwarz : float 364 Schwarz criterion 365 predy_e : array 366 predicted values from reduced form 367 e_pred : array 368 prediction errors using reduced form predicted values 369 pr2 : float 370 Pseudo R squared (squared correlation between y and ypred) 371 pr2_e : float 372 Pseudo R squared (squared correlation between y and ypred_e 373 (using reduced form)) 374 utu : float 375 Sum of squared residuals 376 std_err : array 377 1xk array of standard errors of the betas 378 z_stat : list of tuples 379 z statistic; each tuple contains the pair (statistic, 380 p-value), where each is a float 381 name_y : string 382 Name of dependent variable for use in output 383 name_x : list of strings 384 Names of independent variables for use in output 385 name_w : string 386 Name of weights matrix for use in output 387 name_ds : string 388 Name of dataset for use in output 389 title : string 390 Name of the regression method used 391 392 Examples 393 -------- 394 >>> import numpy as np 395 >>> import libpysal 396 >>> from libpysal.examples import load_example 397 >>> from libpysal.weights import Queen 398 >>> from spreg import ML_Error_Regimes 399 >>> import geopandas as gpd 400 >>> from spreg import ML_Lag 401 >>> np.set_printoptions(suppress=True) #prevent scientific format 402 >>> baltimore = load_example('Baltimore') 403 >>> db = libpysal.io.open(baltimore.get_path("baltim.dbf"),'r') 404 >>> df = gpd.read_file(baltimore.get_path("baltim.shp")) 405 >>> ds_name = "baltim.dbf" 406 >>> y_name = "PRICE" 407 >>> y = np.array(db.by_col(y_name)).T 408 >>> y.shape = (len(y),1) 409 >>> x_names = ["NROOM","NBATH","PATIO","FIREPL","AC","GAR","AGE","LOTSZ","SQFT"] 410 >>> x = np.array([db.by_col(var) for var in x_names]).T 411 >>> w = Queen.from_dataframe(df) 412 >>> w_name = "baltim_q.gal" 413 >>> w.transform = 'r' 414 >>> mllag = ML_Lag(y,x,w,name_y=y_name,name_x=x_names,\ 415 name_w=w_name,name_ds=ds_name) #doctest: +SKIP 416 >>> np.around(mllag.betas, decimals=4) #doctest: +SKIP 417 array([[ 4.3675], 418 [ 0.7502], 419 [ 5.6116], 420 [ 7.0497], 421 [ 7.7246], 422 [ 6.1231], 423 [ 4.6375], 424 [-0.1107], 425 [ 0.0679], 426 [ 0.0794], 427 [ 0.4259]]) 428 >>> "{0:.6f}".format(mllag.rho) #doctest: +SKIP 429 '0.425885' 430 >>> "{0:.6f}".format(mllag.mean_y) #doctest: +SKIP 431 '44.307180' 432 >>> "{0:.6f}".format(mllag.std_y) #doctest: +SKIP 433 '23.606077' 434 >>> np.around(np.diag(mllag.vm1), decimals=4) #doctest: +SKIP 435 array([ 23.8716, 1.1222, 3.0593, 7.3416, 5.6695, 5.4698, 436 2.8684, 0.0026, 0.0002, 0.0266, 0.0032, 220.1292]) 437 >>> np.around(np.diag(mllag.vm), decimals=4) #doctest: +SKIP 438 array([ 23.8716, 1.1222, 3.0593, 7.3416, 5.6695, 5.4698, 439 2.8684, 0.0026, 0.0002, 0.0266, 0.0032]) 440 >>> "{0:.6f}".format(mllag.sig2) #doctest: +SKIP 441 '151.458698' 442 >>> "{0:.6f}".format(mllag.logll) #doctest: +SKIP 443 '-832.937174' 444 >>> "{0:.6f}".format(mllag.aic) #doctest: +SKIP 445 '1687.874348' 446 >>> "{0:.6f}".format(mllag.schwarz) #doctest: +SKIP 447 '1724.744787' 448 >>> "{0:.6f}".format(mllag.pr2) #doctest: +SKIP 449 '0.727081' 450 >>> "{0:.4f}".format(mllag.pr2_e) #doctest: +SKIP 451 '0.7062' 452 >>> "{0:.4f}".format(mllag.utu) #doctest: +SKIP 453 '31957.7853' 454 >>> np.around(mllag.std_err, decimals=4) #doctest: +SKIP 455 array([ 4.8859, 1.0593, 1.7491, 2.7095, 2.3811, 2.3388, 1.6936, 456 0.0508, 0.0146, 0.1631, 0.057 ]) 457 >>> np.around(mllag.z_stat, decimals=4) #doctest: +SKIP 458 array([[ 0.8939, 0.3714], 459 [ 0.7082, 0.4788], 460 [ 3.2083, 0.0013], 461 [ 2.6018, 0.0093], 462 [ 3.2442, 0.0012], 463 [ 2.6181, 0.0088], 464 [ 2.7382, 0.0062], 465 [-2.178 , 0.0294], 466 [ 4.6487, 0. ], 467 [ 0.4866, 0.6266], 468 [ 7.4775, 0. ]]) 469 >>> mllag.name_y #doctest: +SKIP 470 'PRICE' 471 >>> mllag.name_x #doctest: +SKIP 472 ['CONSTANT', 'NROOM', 'NBATH', 'PATIO', 'FIREPL', 'AC', 'GAR', 'AGE', 'LOTSZ', 'SQFT', 'W_PRICE'] 473 >>> mllag.name_w #doctest: +SKIP 474 'baltim_q.gal' 475 >>> mllag.name_ds #doctest: +SKIP 476 'baltim.dbf' 477 >>> mllag.title #doctest: +SKIP 478 'MAXIMUM LIKELIHOOD SPATIAL LAG (METHOD = FULL)' 479 >>> mllag = ML_Lag(y,x,w,method='ord',name_y=y_name,name_x=x_names,\ 480 name_w=w_name,name_ds=ds_name) #doctest: +SKIP 481 >>> np.around(mllag.betas, decimals=4) #doctest: +SKIP 482 array([[ 4.3675], 483 [ 0.7502], 484 [ 5.6116], 485 [ 7.0497], 486 [ 7.7246], 487 [ 6.1231], 488 [ 4.6375], 489 [-0.1107], 490 [ 0.0679], 491 [ 0.0794], 492 [ 0.4259]]) 493 >>> "{0:.6f}".format(mllag.rho) #doctest: +SKIP 494 '0.425885' 495 >>> "{0:.6f}".format(mllag.mean_y) #doctest: +SKIP 496 '44.307180' 497 >>> "{0:.6f}".format(mllag.std_y) #doctest: +SKIP 498 '23.606077' 499 >>> np.around(np.diag(mllag.vm1), decimals=4) #doctest: +SKIP 500 array([ 23.8716, 1.1222, 3.0593, 7.3416, 5.6695, 5.4698, 501 2.8684, 0.0026, 0.0002, 0.0266, 0.0032, 220.1292]) 502 >>> np.around(np.diag(mllag.vm), decimals=4) #doctest: +SKIP 503 array([ 23.8716, 1.1222, 3.0593, 7.3416, 5.6695, 5.4698, 504 2.8684, 0.0026, 0.0002, 0.0266, 0.0032]) 505 >>> "{0:.6f}".format(mllag.sig2) #doctest: +SKIP 506 '151.458698' 507 >>> "{0:.6f}".format(mllag.logll) #doctest: +SKIP 508 '-832.937174' 509 >>> "{0:.6f}".format(mllag.aic) #doctest: +SKIP 510 '1687.874348' 511 >>> "{0:.6f}".format(mllag.schwarz) #doctest: +SKIP 512 '1724.744787' 513 >>> "{0:.6f}".format(mllag.pr2) #doctest: +SKIP 514 '0.727081' 515 >>> "{0:.6f}".format(mllag.pr2_e) #doctest: +SKIP 516 '0.706198' 517 >>> "{0:.4f}".format(mllag.utu) #doctest: +SKIP 518 '31957.7853' 519 >>> np.around(mllag.std_err, decimals=4) #doctest: +SKIP 520 array([ 4.8859, 1.0593, 1.7491, 2.7095, 2.3811, 2.3388, 1.6936, 521 0.0508, 0.0146, 0.1631, 0.057 ]) 522 >>> np.around(mllag.z_stat, decimals=4) #doctest: +SKIP 523 array([[ 0.8939, 0.3714], 524 [ 0.7082, 0.4788], 525 [ 3.2083, 0.0013], 526 [ 2.6018, 0.0093], 527 [ 3.2442, 0.0012], 528 [ 2.6181, 0.0088], 529 [ 2.7382, 0.0062], 530 [-2.178 , 0.0294], 531 [ 4.6487, 0. ], 532 [ 0.4866, 0.6266], 533 [ 7.4775, 0. ]]) 534 >>> mllag.name_y #doctest: +SKIP 535 'PRICE' 536 >>> mllag.name_x #doctest: +SKIP 537 ['CONSTANT', 'NROOM', 'NBATH', 'PATIO', 'FIREPL', 'AC', 'GAR', 'AGE', 'LOTSZ', 'SQFT', 'W_PRICE'] 538 >>> mllag.name_w #doctest: +SKIP 539 'baltim_q.gal' 540 >>> mllag.name_ds #doctest: +SKIP 541 'baltim.dbf' 542 >>> mllag.title #doctest: +SKIP 543 'MAXIMUM LIKELIHOOD SPATIAL LAG (METHOD = ORD)' 544 545 546 """ 547 548 def __init__(self, y, x, w, method='full', epsilon=0.0000001, 549 vm=False, name_y=None, name_x=None, 550 name_w=None, name_ds=None): 551 n = USER.check_arrays(y, x) 552 y = USER.check_y(y, n) 553 USER.check_weights(w, y, w_required=True) 554 x_constant,name_x,warn = USER.check_constant(x,name_x) 555 set_warn(self, warn) 556 method = method.upper() 557 BaseML_Lag.__init__( 558 self, y=y, x=x_constant, w=w, method=method, epsilon=epsilon) 559 # increase by 1 to have correct aic and sc, include rho in count 560 self.k += 1 561 self.title = "MAXIMUM LIKELIHOOD SPATIAL LAG" + \ 562 " (METHOD = " + method + ")" 563 self.name_ds = USER.set_name_ds(name_ds) 564 self.name_y = USER.set_name_y(name_y) 565 self.name_x = USER.set_name_x(name_x, x) 566 name_ylag = USER.set_name_yend_sp(self.name_y) 567 self.name_x.append(name_ylag) # rho changed to last position 568 self.name_w = USER.set_name_w(name_w, w) 569 self.aic = DIAG.akaike(reg=self) 570 self.schwarz = DIAG.schwarz(reg=self) 571 SUMMARY.ML_Lag(reg=self, w=w, vm=vm, spat_diag=False) 572 573def lag_c_loglik(rho, n, e0, e1, W): 574 # concentrated log-lik for lag model, no constants, brute force 575 er = e0 - rho * e1 576 sig2 = spdot(er.T, er) / n 577 nlsig2 = (n / 2.0) * np.log(sig2) 578 a = -rho * W 579 spfill_diagonal(a, 1.0) 580 jacob = np.log(np.linalg.det(a)) 581 # this is the negative of the concentrated log lik for minimization 582 clik = nlsig2 - jacob 583 return clik 584 585def lag_c_loglik_sp(rho, n, e0, e1, I, Wsp): 586 # concentrated log-lik for lag model, sparse algebra 587 if isinstance(rho, np.ndarray): 588 if rho.shape == (1,1): 589 rho = rho[0][0] #why does the interior value change? 590 er = e0 - rho * e1 591 sig2 = spdot(er.T, er) / n 592 nlsig2 = (n / 2.0) * np.log(sig2) 593 a = I - rho * Wsp 594 LU = SuperLU(a.tocsc()) 595 jacob = np.sum(np.log(np.abs(LU.U.diagonal()))) 596 clike = nlsig2 - jacob 597 return clike 598 599def lag_c_loglik_ord(rho, n, e0, e1, evals): 600 # concentrated log-lik for lag model, no constants, Ord eigenvalue method 601 er = e0 - rho * e1 602 sig2 = spdot(er.T, er) / n 603 nlsig2 = (n / 2.0) * np.log(sig2) 604 revals = rho * evals 605 jacob = np.log(1 - revals).sum() 606 if isinstance(jacob, complex): 607 jacob = jacob.real 608 # this is the negative of the concentrated log lik for minimization 609 clik = nlsig2 - jacob 610 return clik 611 612def _test(): 613 import doctest 614 start_suppress = np.get_printoptions()['suppress'] 615 np.set_printoptions(suppress=True) 616 doctest.testmod() 617 np.set_printoptions(suppress=start_suppress) 618 619