1""" 2Diagnostics for SUR and 3SLS estimation 3""" 4 5__author__= "Luc Anselin lanselin@gmail.com, \ 6 Pedro V. Amaral pedrovma@gmail.com \ 7 Tony Aburaad taburaad@uchicago.edu" 8 9 10import numpy as np 11import scipy.stats as stats 12import numpy.linalg as la 13from .sur_utils import sur_dict2mat,sur_mat2dict,sur_corr,spdot 14from .regimes import buildR1var,wald_test 15 16 17__all__ = ['sur_setp','sur_lrtest','sur_lmtest','lam_setp','surLMe','surLMlag'] 18 19 20def sur_setp(bigB,varb): 21 ''' 22 Utility to compute standard error, t and p-value 23 24 Parameters 25 ---------- 26 bigB : dictionary 27 of regression coefficient estimates, 28 one vector by equation 29 varb : array 30 variance-covariance matrix of coefficients 31 32 Returns 33 ------- 34 surinfdict : dictionary 35 with standard error, t-value, and 36 p-value array, one for each equation 37 38 ''' 39 vvb = varb.diagonal() 40 n_eq = len(bigB.keys()) 41 bigK = np.zeros((n_eq,1),dtype=np.int_) 42 for r in range(n_eq): 43 bigK[r] = bigB[r].shape[0] 44 b = sur_dict2mat(bigB) 45 se = np.sqrt(vvb) 46 se.resize(len(se),1) 47 t = np.divide(b,se) 48 tp = stats.norm.sf(abs(t))*2 49 surinf = np.hstack((se,t,tp)) 50 surinfdict = sur_mat2dict(surinf,bigK) 51 return surinfdict 52 53def lam_setp(lam,vm): 54 """ 55 Standard errors, t-test and p-value for lambda in SUR Error ML 56 57 Parameters 58 ---------- 59 lam : array 60 n_eq x 1 array with ML estimates for spatial error 61 autoregressive coefficient 62 vm : array 63 n_eq x n_eq subset of variance-covariance matrix for 64 lambda and Sigma in SUR Error ML 65 (needs to be subset from full vm) 66 67 Returns 68 ------- 69 : tuple 70 with arrays for standard error, t-value and p-value 71 (each element in the tuple is an n_eq x 1 array) 72 73 """ 74 vvb = vm.diagonal() 75 se = np.sqrt(vvb) 76 se.resize(len(se),1) 77 t = np.divide(lam,se) 78 tp = stats.norm.sf(abs(t))*2 79 return (se,t,tp) 80 81def sur_lrtest(n,n_eq,ldetS0,ldetS1): 82 ''' 83 Likelihood Ratio test on off-diagonal elements of Sigma 84 85 Parameters 86 ---------- 87 n : int 88 cross-sectional dimension (number of observations for an equation) 89 n_eq : int 90 number of equations 91 ldetS0 : float 92 log determinant of Sigma for OLS case 93 ldetS1 : float 94 log determinant of Sigma for SUR case (should be iterated) 95 96 Returns 97 ------- 98 (lrtest,M,pvalue) : tuple 99 with value of test statistic (lrtest), 100 degrees of freedom (M, as an integer) 101 p-value 102 103 ''' 104 M = n_eq * (n_eq - 1)/2.0 105 lrtest = n * (ldetS0 - ldetS1) 106 pvalue = stats.chi2.sf(lrtest,M) 107 return (lrtest,int(M),pvalue) 108 109 110def sur_lmtest(n,n_eq,sig): 111 ''' 112 Lagrange Multiplier test on off-diagonal elements of Sigma 113 114 Parameters 115 ---------- 116 n : int 117 cross-sectional dimension (number of observations for an equation) 118 n_eq : int 119 number of equations 120 sig : array 121 inter-equation covariance matrix for null model (OLS) 122 123 Returns 124 ------- 125 (lmtest,M,pvalue) : tuple 126 with value of test statistic (lmtest), 127 degrees of freedom (M, as an integer) 128 p-value 129 ''' 130 R = sur_corr(sig) 131 tr = np.trace(np.dot(R.T,R)) 132 M = n_eq * (n_eq - 1)/2.0 133 lmtest = (n/2.0) * (tr - n_eq) 134 pvalue = stats.chi2.sf(lmtest,M) 135 return (lmtest,int(M),pvalue) 136 137 138def surLMe(n_eq,WS,bigE,sig): 139 """ 140 Lagrange Multiplier test on error spatial autocorrelation in SUR 141 142 Parameters 143 ---------- 144 n_eq : int 145 number of equations 146 WS : array 147 spatial weights matrix in sparse form 148 bigE : array 149 n x n_eq matrix of residuals by equation 150 sig : array 151 cross-equation error covariance matrix 152 153 Returns 154 ------- 155 (LMe,n_eq,pvalue) : tuple 156 with value of statistic (LMe), degrees 157 of freedom (n_eq) and p-value 158 159 """ 160 # spatially lagged residuals 161 WbigE = WS * bigE 162 # score 163 EWE = np.dot(bigE.T,WbigE) 164 sigi = la.inv(sig) 165 SEWE = sigi * EWE 166 #score = SEWE.sum(axis=1) 167 #score.resize(n_eq,1) 168 # note score is column sum of Sig_i * E'WE, a 1 by n_eq row vector 169 # previously stored as column 170 score = SEWE.sum(axis=0) 171 score.resize(1,n_eq) 172 173 174 # trace terms 175 WW = WS * WS 176 trWW = np.sum(WW.diagonal()) 177 WTW = WS.T * WS 178 trWtW = np.sum(WTW.diagonal()) 179 # denominator 180 SiS = sigi * sig 181 Tii = trWW * np.identity(n_eq) 182 tSiS = trWtW * SiS 183 denom = Tii + tSiS 184 idenom = la.inv(denom) 185 # test statistic 186 #LMe = np.dot(np.dot(score.T,idenom),score)[0][0] 187 # score is now row vector 188 LMe = np.dot(np.dot(score,idenom),score.T)[0][0] 189 pvalue = stats.chi2.sf(LMe,n_eq) 190 return (LMe,n_eq,pvalue) 191 192def surLMlag(n_eq,WS,bigy,bigX,bigE,bigYP,sig,varb): 193 """ 194 Lagrange Multiplier test on lag spatial autocorrelation in SUR 195 196 Parameters 197 ---------- 198 n_eq : int 199 number of equations 200 WS : spatial weights matrix in sparse form 201 bigy : dictionary 202 with y values 203 bigX : dictionary 204 with X values 205 bigE : array 206 n x n_eq matrix of residuals by equation 207 bigYP : array 208 n x n_eq matrix of predicted values by equation 209 sig : array 210 cross-equation error covariance matrix 211 varb : array 212 variance-covariance matrix for b coefficients (inverse of Ibb) 213 214 Returns 215 ------- 216 (LMlag,n_eq,pvalue) : tuple 217 with value of statistic (LMlag), degrees 218 of freedom (n_eq) and p-value 219 220 """ 221 # Score 222 Y = np.hstack((bigy[r]) for r in range(n_eq)) 223 WY = WS * Y 224 EWY = np.dot(bigE.T,WY) 225 sigi = la.inv(sig) 226 SEWE = sigi * EWY 227 score = SEWE.sum(axis=0) # column sums 228 score.resize(1,n_eq) # score as a row vector 229 230 # I(rho,rho) as partitioned inverse, eq 72 231 # trace terms 232 WW = WS * WS 233 trWW = np.sum(WW.diagonal()) #T1 234 WTW = WS.T * WS 235 trWtW = np.sum(WTW.diagonal()) #T2 236 237 # I(rho,rho) 238 SiS = sigi * sig 239 Tii = trWW * np.identity(n_eq) #T1It 240 tSiS = trWtW * SiS 241 firstHalf = Tii + tSiS 242 WbigYP = WS * bigYP 243 inner = np.dot(WbigYP.T, WbigYP) 244 secondHalf= sigi * inner 245 Ipp= firstHalf + secondHalf #eq. 75 246 247 # I(b,b) inverse is varb 248 249 # I(b,rho) 250 bp = sigi[0,] * spdot(bigX[0].T,WbigYP) #initialize 251 for r in range(1,n_eq): 252 bpwork = sigi[r,] * spdot(bigX[r].T,WbigYP) 253 bp = np.vstack((bp,bpwork)) 254 # partitioned part 255 i_inner = Ipp - np.dot(np.dot(bp.T, varb), bp) 256 # partitioned inverse of information matrix 257 Ippi = la.inv(i_inner) 258 259 # test statistic 260 LMlag = np.dot(np.dot(score,Ippi),score.T)[0][0] 261 # p-value 262 pvalue = stats.chi2.sf(LMlag,n_eq) 263 return (LMlag,n_eq,pvalue) 264 265def sur_chow(n_eq,bigK,bSUR,varb): 266 """ 267 test on constancy of regression coefficients across equations in 268 a SUR specification 269 270 Note: requires a previous check on constancy of number of coefficients 271 across equations; no other checks are carried out, so it is possible 272 that the results are meaningless if the variables are not listed in 273 the same order in each equation. 274 275 Parameters 276 ---------- 277 n_eq : int 278 number of equations 279 bigK : array 280 with the number of variables by equation (includes constant) 281 bSUR : dictionary 282 with the SUR regression coefficients by equation 283 varb : array 284 the variance-covariance matrix for the SUR regression 285 coefficients 286 287 Returns 288 ------- 289 test : array 290 a list with for each coefficient (in order) a tuple with the 291 value of the test statistic, the degrees of freedom, and the 292 p-value 293 294 """ 295 kr = bigK[0][0] 296 test = [] 297 bb = sur_dict2mat(bSUR) 298 kf = 0 299 nr = n_eq 300 df = n_eq - 1 301 for i in range(kr): 302 Ri = buildR1var(i,kr,kf,0,nr) 303 tt,p = wald_test(bb,Ri,np.zeros((df,1)),varb) 304 test.append((tt,df,p)) 305 return test 306 307def sur_joinrho(n_eq,bigK,bSUR,varb): 308 """ 309 Test on joint significance of spatial autoregressive coefficient in SUR 310 311 Parameters 312 ---------- 313 n_eq : int 314 number of equations 315 bigK : array 316 n_eq x 1 array with number of variables by equation 317 (includes constant term, exogenous and endogeneous and 318 spatial lag) 319 bSUR : dictionary 320 with regression coefficients by equation, with 321 the spatial autoregressive term as last 322 varb : array 323 variance-covariance matrix for regression coefficients 324 325 Returns 326 ------- 327 : tuple 328 with test statistic, degrees of freedom, p-value 329 330 """ 331 bb = sur_dict2mat(bSUR) 332 R = np.zeros((n_eq,varb.shape[0])) 333 q = np.zeros((n_eq,1)) 334 kc = -1 335 for i in range(n_eq): 336 kc = kc + bigK[i] 337 R[i,kc] = 1 338 w,p = wald_test(bb,R,q,varb) 339 return (w,n_eq,p) 340